{"cells":[{"metadata":{},"cell_type":"markdown","source":"## <u>About this kernel</u>\n\nIn this kernel, I propose following two methods. First, I will explain outline of the methods.\n\n**1. For calculate Emissions Factor**\n\n- I use NO2_column_number_density of Sentinel-5P dataset.\n- I regard one data (such as s5p_no2_20180701T161259_20180707T175356.tif) as one day's worth. If na is contained, fill na by mean of the data.\n- For calculate amount of NOx, sum up the density in focused area.\n- For calculate total estimated generation in focused area, sum up estimated generation from \"Global Power Plant Database\" in the area.\n- For calculate emissions factor, devide the total NOx amount by total estimated generation.\n\n**2. For calculate Merginal Emissions Factor**\n\n- Create dataframe of NO2_column_number_density, longtitude, latitude, wind propaty.\n- By k-means divide data of NO2_column_number_density into some regional groups using the dataframe.\n- Using method 1 and some additional assumptions, I calculate Emissions Factors of each power plants. \n  Then, sum the NOx amount and the estimated power generation only within the group to which the power plant belongs.\n- Using the Emissions Factors calculate merginal emissions factor.\n\nI'll explain the details of the method using NO2 as the theme."},{"metadata":{},"cell_type":"markdown","source":"## Why pay attention to NO2?\n\nThe simplest way to calculate emission factors is dividing the actual measured value of the amount of NO2 emitted by the power plant by the amount of power generated per hour. This value seems be called \"actual emissions factor\". (I refered P.30 of [1]　for this actual emissions factor.)\n\nBut since NOx reacts with chemicals in the atmosphere and changes into various forms, the amount measured at the time of generation is not always the same as the amount of NO2 actually present in the environment. \n\nTherefore, if the NO2 emissions factor can be calculated accurately, it will be a reference for other gases."},{"metadata":{},"cell_type":"markdown","source":"## About calculation\n\nI aimed to calculate g/GWh unit Emissions Factor (such as g/GWh, T/GWh). In some document I read, this unit seems standard.\n\nI used fllowing fomula.\n\n### Emissions(mol) = NO2_column_number_density(mol/m^2) \\* area(m^2)\n### Emissons(T) = Emissions(mol) \\* 48(g/mol) \\* 1e^-6\n### Emissons Factor(T/GWh) = Emissons(T) / Estimated Generation(GWh)"},{"metadata":{},"cell_type":"markdown","source":"## Load libraries"},{"metadata":{"trusted":false},"cell_type":"code","source":"from decimal import Decimal, ROUND_HALF_UP, ROUND_HALF_EVEN\nfrom datetime import datetime, timedelta\nimport folium\nimport glob\nimport math\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport os\nimport pandas as pd\nimport rasterio as rio\nimport seaborn as sns\nfrom sklearn.cluster import KMeans\nimport tifffile as tiff ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Utilities\n\nI put it all together tools that are not related to the subject but are required."},{"metadata":{"trusted":false},"cell_type":"code","source":"#If you want to analyze other area, please change here!\nLAT_MAX = 18.563112930177304\nLAT_MIN = 17.903121359128956\nLON_MAX = -65.19437297127342\nLON_MIN = -67.32297404549217","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"When trying other places, you need to change the following constants: The constants represent the latitude and longitude range of the data."},{"metadata":{"_kg_hide-input":true,"trusted":false},"cell_type":"code","source":"def overlay_image_on_puerto_rico_df(df, img, zoom):\n    \"\"\"\n    show image on google map with marker of power plants.\n    \"\"\"\n    lat_map=df.iloc[[0]].loc[:,[\"latitude\"]].iat[0,0]\n    lon_map=df.iloc[[0]].loc[:,[\"longitude\"]].iat[0,0]\n    m = folium.Map([lat_map, lon_map], zoom_start=zoom, tiles= 'Stamen Terrain')\n    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\n    folium.raster_layers.ImageOverlay(\n        image=img,\n        bounds = [[LAT_MAX,LON_MIN,],[LAT_MIN,LON_MAX]],\n        #bounds = [[18.56,-67.32,],[17.90,-65.194]],\n        colormap=lambda x: (1, 0, 0, x),\n    ).add_to(m)\n    \n    for i in range(0,len(df)):\n        popup = folium.Popup(str(df.primary_fuel[i:i+1]))\n        folium.Marker([df[\"latitude\"].iloc[i],df[\"longitude\"].iloc[i]],\n                     icon=folium.Icon(icon_color='red',icon ='bolt',prefix='fa',color=color[df.primary_fuel.iloc[i]])).add_to(m)\n        \n    return m","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":false},"cell_type":"code","source":"def split_column_into_new_columns(dataframe,column_to_split,new_column_one,begin_column_one,end_column_one):\n    \"\"\"\n    Add latitude and longtitude to dataframe.\n    \"\"\"\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":{},"cell_type":"markdown","source":"# Data overview\n\nFor the first point, we get power plants data from \"Global Power Plant Database\" in the area.\n\nIf we calculate emissions factor in other place, we have to get data of power plants in coressponding area and replace here."},{"metadata":{"trusted":false},"cell_type":"code","source":"#import Global Power Plant Database data\npower_plants = pd.read_csv('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\n\n#Make latitude and longitude data easier to use\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() \n\n#sort data by their capacity\npower_plants_df = power_plants.sort_values('capacity_mw',ascending=False).reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"power_plants_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are 35 power plants in Puerto Rico.\n\nHydrogen, solar and wind are NOx free power source. \n\nOil, coal, and gas emits NOx. \n\nAbout 45% power plants in Puerto Rico emit NOx."},{"metadata":{"trusted":false},"cell_type":"code","source":"#Check NO2_column_number_density of Sentinel-5P Data\nimage = tiff.imread('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif')\noverlay_image_on_puerto_rico_df(power_plants_df,image[:,:,0],8)\n\n##From https://www.kaggle.com/paultimothymooney/explore-image-metadata-s5p-gfs-gldas\n##You can check which bands correspond which properties. \n##Now, band1 is NO2_column_number_density.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Correspondence of marker color and fuel are below.\n\n- 'Hydrogen' : 'lightblue' \n- 'Solar' : 'orange'\n- 'Oil' : 'darkblue'\n- 'Coal' : 'black'\n- 'Gas' : 'lightgray' \n- 'Wind' : 'green'\n\nOvarlayed data are (148, 475) shape numpy arrays. Area of 1 pixel is 0.25 m^2※.\n\n※Use following values to lead this value. The original data is located at latitudes 17.9 to 18.6 degrees and longitudes -67.3 to -65.2 degrees. And longitude 1degree is 111km."},{"metadata":{},"cell_type":"markdown","source":"# 1. For calculate Emissions Factor\n\nNow, I try to calculate emissions factor whole Puerto Rico.\n\nFirst, I calculate total estimated generation. Import and pre-process NO2_column_number_density data."},{"metadata":{},"cell_type":"markdown","source":"## Preparation "},{"metadata":{"trusted":false},"cell_type":"code","source":"#Calculate total estimated electricity generation(GWh)\nquantity_of_electricity_generated = np.sum(power_plants_df['estimated_generation_gwh'])\nprint('Quanity of Electricity Generated: ', quantity_of_electricity_generated)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"#import path of Sentinel-5P Data\ns5p_no2_timeseries = glob.glob('../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/*')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"NO2_column_number_density is band1 of Sentinel-5P dataset. If you want to calculate emissions factor in other place, please take data of the region from Sentinel-5P dataset and replace it at here.\n\nIf you check sample s5p_no2_timeseries, you can find that there are data which starts same date. And they effect result (values become abnormally big because of duplication). So I use most early data in each days simply."},{"metadata":{"trusted":false},"cell_type":"code","source":"s5p_no2_timeseries_no_duplication = []\nchecked_date = []\n\nfor data in sorted(s5p_no2_timeseries):\n     \n    data_date =  datetime.strptime(data[:79], '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_%Y%m%d')\n    data_date = data_date.strftime(\"%Y/%m%d\")\n    \n    if not data_date in checked_date:\n        checked_date.append(data_date)\n        s5p_no2_timeseries_no_duplication.append(data)\n\n#Data path without duplicates\ns5p_no2_timeseries = s5p_no2_timeseries_no_duplication","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Caluculate monthly emissions factor\n\nSecond, I calculate emissions factor monthly. By this way, we can calculate emissions factor in monthly resolution."},{"metadata":{"trusted":false},"cell_type":"code","source":"#Divide the data by month.\ndata_monthly_divided = {}\nfor data in s5p_no2_timeseries:\n     \n    data_date =  datetime.strptime(data[:77], '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_%Y%m')\n    data_date = data_date.strftime(\"%Y/%m\")\n    \n    if not data_date in data_monthly_divided.keys():\n        data_monthly_divided[data_date] = []\n        \n    data_monthly_divided[data_date].append(data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"month = []\nemissions = []\n\nfor key in sorted(data_monthly_divided.keys()):\n    total_emissions = []\n    datas = data_monthly_divided[key]\n    \n    for data in datas:        \n        img = tiff.imread(data)[:,:,0] #import data here\n        img = np.nan_to_num(img, nan=np.nanmean(img))  #fill nan by average  \n        total_emission = np.nansum(img)  #take total NO2 density of the data\n        \n        total_emissions.append(total_emission)\n    \n    #take monthly total density.\n    month.append(key)\n    emissions.append(np.nansum(total_emissions))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We import data here.\n\n> img = tiff.imread(data)[:,:,0] #import data here\n\nJust change imported data, we use other place and scale data. "},{"metadata":{"trusted":false},"cell_type":"code","source":"#calculate amount of NO2.\n#amount[T] = density[mol/m^2] * 0.25m^2 * number of whole pixels * 46.0055[g/mol] * 1e-6\n\nemissions = np.array(emissions) * ((0.25 * img.shape[0]*img.shape[1]) * 46.0055 *1e-6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"results_monthly = pd.DataFrame(columns=['month', 'emission','emisson factor'])\nresults_monthly = pd.DataFrame({'emission':emissions,\n                       'emission factor':emissions/(quantity_of_electricity_generated)}, #devide emissions by estimated generation\n                    index=month)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"results_monthly.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"fig = plt.figure(figsize=(30, 4))\nax = results_monthly[\"emission factor\"].plot()\nplt.title('Monthly Emissions Factor in Puerto Rico')\nax.set(xlabel='YYYY/mm', ylabel='Emission factor [T/GWh]')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Caluculate minimal time span (Appendix)\n\nSentinel-5P starts measuring every day, so using Sentinel-5P data, we can calculate more detailed values.\n\nThird, I try to calculate emissions factor in this minimal time span for appendix."},{"metadata":{"trusted":false},"cell_type":"code","source":"#Divide the data by minimal time span.\ndata_minimal_divided = {}\nfor data in s5p_no2_timeseries:\n     \n    data_date =  datetime.strptime(data[:79], '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_%Y%m%d')\n    data_date = data_date.strftime(\"%Y/%m/%d\")\n    \n    if not data_date in data_minimal_divided.keys():\n        data_minimal_divided[data_date] = [data]","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"span = []\nemissions = []\n\nfor key in sorted(data_minimal_divided.keys()):\n    total_emissions = []\n    datas = data_minimal_divided[key]\n    for data in datas:\n        \n        img = tiff.imread(data)[:,:,0]\n        img = np.nan_to_num(img, nan=np.nanmean(img))\n        total_emission = np.nansum(img)\n        \n        total_emissions.append(total_emission)\n    \n    span.append(key)\n    emissions.append(np.nansum(total_emissions))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"#calculate amount of NO2.\n#amount[T] = density[mol/m^2] * 0.25m^2 * number of whole pixels * 46.0055[g/mol] * 1e-6\nemissions = np.array(emissions) * ((0.25 * img.shape[0]*img.shape[1]) * 46.0055 *1e-6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"results_minimal = pd.DataFrame(columns=['week', 'emission','emisson factor'])\nresults_minimal = pd.DataFrame({'emission':emissions,\n                       'emission factor':emissions/(quantity_of_electricity_generated)},\n                    index=span)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"results_minimal.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"fig = plt.figure(figsize=(30, 4))\nax = results_minimal[\"emission factor\"].plot()\nplt.title('Monthly Mean Simplified Emissions Factor in Puerto Rico')\nax.set(xlabel='YYYY/mm/dd', ylabel='Emission factor [T/GWh]')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Are these estimates correct?"},{"metadata":{"trusted":false},"cell_type":"code","source":"print( \"Total NO2 emissions in Puerto Rico is\", int(np.sum(emissions) * 12), \"T/year\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I think this value is not far away from real value. According to page21 of [2], total NO2 emissions from factories and business sites of Chiba prefecture in Japan is 41944 T/year. Puerto Rico is 2.5 times larger than Chiba, but Chiba has very big industrial zone."},{"metadata":{},"cell_type":"markdown","source":"# 2. For calculate Merginal Emissions Factor"},{"metadata":{},"cell_type":"markdown","source":"## Overview of objectives and methods\nThere are 35 power plants and 6 fuel types throughout Puerto Rico. The capacity for each fuel type can be summarized in the graph below."},{"metadata":{"_kg_hide-input":true,"trusted":false},"cell_type":"code","source":"#From https://www.kaggle.com/ajulian/capacity-factor-in-power-plants\n\ntotal_capacity_mw = power_plants_df['capacity_mw'].sum()\nprint('Total Installed Capacity: '+'{:.2f}'.format(total_capacity_mw) + ' MW')\ncapacity = (power_plants_df.groupby(['primary_fuel'])['capacity_mw'].sum()).to_frame()\ncapacity = capacity.sort_values('capacity_mw',ascending=False)\ncapacity['percentage_of_total'] = (capacity['capacity_mw']/total_capacity_mw)*100\nax = capacity.sort_values(by='percentage_of_total', ascending=True)['percentage_of_total'].plot(kind='bar',color=['lightblue', 'green', 'orange', 'black','lightgray','darkblue'])\nax.set(ylabel='percentage')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since the marginal emissions factor is an emissions factor of the power plant that supplies power at a certain time, it is necessary to estimate the NOx emissions for each power plant. The sample data is provided for the entire Puerto Rico, but since Puerto Rico is large, it is not reasonable to simply distribute the total emissions by power generation. In addition, Google Earth Engine can collect data in a fixed form such as a rectangle, but the distribution of NOx emitted from a certain power plant changes dynamically depending on the characteristics such as land, weather and so on. Therefore, there is a need for a more flexible attribution method of which NOx is emitted from which power plant. \n\nI tought that it is nice idea to divide whole data into some geographical groups by clustering method such as k-means. \n\nIdeally I want to devide whole data into groups as many as power plants in Puerto Rico. But it is difficult because there are some power plants very near each other. \n\nWhen there are some power plants in each group, I devide emissions by multiplicatting capacity_mw rate※ in the area. I think capacity reflect power plants's scale of NOx emissions.\n\n※capacity_mw rate = (capacity_mw of specific plant)/  (capacity_mw of all plants in the area)"},{"metadata":{},"cell_type":"markdown","source":"## Can we devide only NO2 density information?\n\nFirst, I confirm NO2 distribution, and whether we should devide the data by which only NO2 density but also other propaties.\n\nI roughly classified NO2 distribution by k-means only using its own data."},{"metadata":{"trusted":false},"cell_type":"code","source":"image = tiff.imread('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif')\nmonotonous = KMeans(n_clusters=3, random_state=6).fit_predict(image[:,:,0].reshape(-1, 1))\noverlay_image_on_puerto_rico_df(power_plants_df, monotonous.reshape((148, 475)), 8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I classified the data into three levels so that the darker the color, the higher the NOx concentration.\n\nAs you can see from this figure, the NO2 gas concentration is not easy to understand distribution.\n\nThe gas is not distributed like a circle or band from a power plant. It seems be mottled.\n\nSo it seems necessary to include not only density but also other effects, in order to divide data into groups robustly and accurately."},{"metadata":{},"cell_type":"markdown","source":"## Divide into groups with wind feature\n\nI tried to devide data into groups including NO2 density and wind feature.\n\nSince NO2 gases are flowed by wind, and accumulate place where no wind, I considered that I should include this propaties to get robust and accurate result.\n\nWe can access following by Global Forecast System 384-Hour Predicted Atmosphere Data,\n\n- u_component_of_wind_10m_above_ground\n\n- v_component_of_wind_10m_above_ground\n\nand by GLDAS-2.1: Global Land Data Assimilation System,\n\n- Wind_f_inst\n\nIf you want to calculate for other region, please replace here."},{"metadata":{"trusted":false},"cell_type":"code","source":"gldas_files = glob.glob('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gldas/*')\ngldas_files = sorted(gldas_files)\ngfs_files = glob.glob('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/*')\ngfs_files = sorted(gfs_files)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For simplicity, I focus on the data from 7/1 to 7/7."},{"metadata":{"trusted":false},"cell_type":"code","source":"#Data are separated by day for the future, but it is not necessary.\ngldas_files_par_day = []\nfor i in range(0,len(gldas_files[6:54]),8):\n    gldas_files_par_day.append(gldas_files[i:i+8])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"#Data are separated by day for the future, but it is not necessary.\ngfs_files_par_day = []\nfor i in range(0,len(gfs_files[3:27]),4):\n    #print(gfs_files[i:i+4])\n    gfs_files_par_day.append(gfs_files[i:i+4])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Create input for clustering model\n\nGLDAS data are taken by 3 hour interval and GFS data are taken by 6 hour interval. So we have to adjust time span to Sentinel-5P OFFL NO2 data."},{"metadata":{"trusted":false},"cell_type":"code","source":"ave_wind_u = []\nave_wind_v = []\nave_wind_speed = []\n\n#Get data of a day\nfor i in range(len(gfs_files_par_day)):\n    gfs_tmp = gfs_files_par_day[i]\n    gldas_tmp = gldas_files_par_day[i]\n    array_wind_u = []\n    array_wind_v = []\n    array_wind_speed = []\n    \n    #Get datas in the day\n    for j in range(len(gfs_tmp)):\n        gfs_image_u = tiff.imread(gfs_tmp[j])[:,:,3]\n        gfs_image_v = tiff.imread(gfs_tmp[j])[:,:,4]\n        gldas_image1 = tiff.imread(gldas_tmp[2*j])[:,:,11]\n        gldas_image2 = tiff.imread(gldas_tmp[2*j + 1])[:,:,11]\n\n        #fill na by mean\n        gfs_image_u = np.nan_to_num(gfs_image_u, nan=np.nanmean(gfs_image_u))\n        gfs_image_v = np.nan_to_num(gfs_image_v, nan=np.nanmean(gfs_image_v))\n        gldas_image1 = np.nan_to_num(gldas_image1, nan=np.nanmean(gldas_image1))\n        gldas_image2 = np.nan_to_num(gldas_image2, nan=np.nanmean(gldas_image2))\n        \n        #GLDAS has twice detailed time span than GFS\n        gldas_image = (gldas_image1 + gldas_image2)/2\n        \n        array_wind_u.append(gfs_image_u)\n        array_wind_v.append(gfs_image_v)\n        array_wind_speed.append(gldas_image)\n    \n#Calculate average        \nave_wind_u = np.nanmean(np.array(array_wind_u), axis=0)\nave_wind_v = np.nanmean(np.array(array_wind_v), axis=0)\nave_wind_speed = np.nanmean(np.array(array_wind_speed), axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"lon = []\nlat = []\nNO2 = []\nwind_u = []\nwind_v = []\nwind_speed = []\n\nfor i in range(image[:,:,0].shape[0]):\n    for j in range(image[:,:,0].shape[1]):\n        #print(image[:,:,0][i,j])\n        NO2.append(image[:,:,0][i,j])\n        lon.append(i)\n        lat.append(j)\n        wind_u.append(ave_wind_u.reshape((148, 475))[i,j])\n        wind_v.append(ave_wind_v.reshape((148, 475))[i,j])\n        wind_speed.append(ave_wind_speed.reshape((148, 475))[i,j])\n        \nNO2 = np.array(NO2)\nlon = np.array(lon)\nlat = np.array(lat)\nwind_u = np.array(wind_u)\nwind_v = np.array(wind_v)\nwind_spped = np.array(wind_speed)\n        \nfeatures_df = pd.DataFrame(columns=['NO2', 'lat', 'lon', 'wind_u', 'wind_v', 'wind_speed'])\nfeatures_df = pd.DataFrame({\n                    'NO2': NO2/max(NO2),\n                    'lat': lat/max(lat),\n                    'lon': lon/max(lon),\n                    'wind_u' : wind_u/(- min(wind_u)),\n                    'wind_v' : wind_v/(- min(wind_v)),\n                    'wind_speed': wind_speed/max(wind_speed)})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Let's divide into groups\n\nI'll divide whole NO2 data to several groups by k-means and created data. By k-means, we have to decide the number of clusters. To decide this, I check the map of Puerto Rico again."},{"metadata":{"trusted":false},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_df, np.zeros((148, 475)), 8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Power plants using oil, coal and gas fuel emits NO2, so I focus of these power plants. From the geographical features and distance perspective, I decide following 7 group.\n  1. Upper left of main island\n  2. Upper right of main island\n  3. Left of main island\n  4. Right of main island\n  5. Lower left of main island\n  6. Lower right of main island\n  7. Vieques island\n\nAs you can see, we have to decide the number of clusters manually. If we use data that is too geographically large, this feature may be a problem.However, there is no problem with the size of subnations."},{"metadata":{"trusted":false},"cell_type":"code","source":"group_pred = KMeans(n_clusters=7, random_state=3).fit_predict(features_df)\nplt.figure()\nsns.heatmap(group_pred.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_df, group_pred.reshape((148, 475)), 8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfortunately, the central mountain area has become one group.😅 However, it seems to reflect the geographic features well. I think that the wind propaties implicitly incorporate geographic features into clustering."},{"metadata":{},"cell_type":"markdown","source":"## Calculate emissions factor of each power plants\n\nNext step, I caclulate that each power plants can generate what percentage of electricity in the group. This is because I want to estimate each power plant's NO2 emissions by multipling the rate to total emissions in the group. Since capacity reflect the scale of each power plant, I measure this rate by capacity. This is just assumption, so we have to devise more good index."},{"metadata":{"trusted":false},"cell_type":"code","source":"def which_pixel_and_group(df,pred,img):\n    \"\"\"\n    Add information (which pixel and class label) to input DataFrame\n    \n    Parameters\n    ----------\n    df : pandas.DataFrame\n        This dataframe must have latitude and longitude.\n    pred : numpy.array\n        Label classified by k-means.\n    img : numpy.array\n        NO2 density.\n\n    Returns\n    -------\n    df : pandas.DataFrame\n        Information of powerplant added pixel and class label.\n    \"\"\"\n    \n    lat_pixel = []\n    lon_pixel = []\n    kmean_groups = []\n    \n    for i in range(len(df)):\n        lat = float(df.iloc[[i]].loc[:,[\"latitude\"]].iat[0,0])\n        lon = float(df.iloc[[i]].loc[:,[\"longitude\"]].iat[0,0])\n\n    \n        f_lat = (lat - LAT_MIN)*img.shape[0]/(LAT_MAX - LAT_MIN)\n        f_lon = (lon + LON_MAX)*img.shape[1]/(-LON_MIN + LON_MAX)\n        f_lat_int = int(Decimal(str(f_lat -1)).quantize(Decimal('0'), rounding=ROUND_HALF_UP))\n        f_lon_int = int(Decimal(str(f_lon -1)).quantize(Decimal('0'), rounding=ROUND_HALF_UP))\n        pixel = (f_lat_int - 1) * img.shape[1] + f_lon_int\n        \n        lat_pixel.append(img.shape[0] - f_lat_int)\n        lon_pixel.append(f_lon_int)\n        kmean_groups.append(pred[pixel])\n \n        \n    df[\"lat_pixel\"] = lat_pixel\n    df[\"lon_pixel\"] = lon_pixel\n    df[\"kmean_group\"] = kmean_groups\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def calc_gwh_in_group(df):\n    \"\"\"\n    Add information (which pixel and class label) to input DataFrame\n    \n    Parameters\n    ----------\n    df : pandas.DataFrame\n        This dataframe must have primary_fuel, kmean_group and estimated_generation_gwh.\n\n    Returns\n    -------\n    df : pandas.DataFrame\n        Information of powerplant added gwh_rate_group.\n    \"\"\"\n    \n    pplant_gwh_rate_ingroup = []\n    \n    for i in range(len(df)):\n        \n        #Exclude power plants that do not emit NO2\n        if not df.iloc[[i]].loc[:,[\"primary_fuel\"]].iat[0,0] in [\"Oil\",\"Gas\", \"Oil\"]:\n            pplant_gwh_rate_ingroup.append(0)\n            \n        else:      \n            pplant_cap = df.iloc[[i]].loc[:,[\"capacity_mw\"]].iat[0,0]\n            pplant_group = df.iloc[[i]].loc[:,[\"kmean_group\"]].iat[0,0]\n        \n            pplants_emitsno2_group = df[(df[\"kmean_group\"]==pplant_group) | \\\n                                         power_plants_df[\"primary_fuel\"].map(lambda primary_fuel: primary_fuel in [\"Oil\",\"Gas\", \"Oil\"])]\n        \n            total_cap_ingroup = sum(pplants_emitsno2_group[\"capacity_mw\"])\n            pplant_gwh_rate_ingroup.append(pplant_cap/total_cap_ingroup)\n        \n    df[\"cap_rate_group\"] = pplant_gwh_rate_ingroup\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"power_plants_df = which_pixel_and_group(power_plants_df, group_pred, image)\npower_plants_df = calc_gwh_in_group(power_plants_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"power_plants_df.loc[:,[\"name\", \"primary_fuel\", \"kmean_group\",\"cap_rate_group\"]].head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Calculate each own emission factor of power plants\n\nLast, I calculate emissions factor of each power plants (own_EF in following code). To calculate the values, I calculate total NO2 amount in each group (no2amount_dict_group), and multiple the values to \"cap_rate_group\"."},{"metadata":{"trusted":false},"cell_type":"code","source":"def calc_no2amount_each_group(img, pred):\n    \"\"\"\n    Calculate amount of substance of NO2 for each group.\n    \n    Parameters\n    ----------\n    img : numpy.array\n        NO2 density.\n    pred : numpy.array\n        Label classified by k-means.\n\n    Returns\n    -------\n    no2amount : dictionary\n        value: group number, value: amount of substance of NO2(g)\n    \"\"\"\n    no2amount = dict()\n    \n    for i in set(pred):\n        no2amount[i] = 0\n    \n    pred = pred.reshape(img[:,:,0].shape)\n    \n    for i in range(img.shape[0]):\n        for j in range(img.shape[1]):\n            no2amount[pred[i,j]] += img[:,:,0][i,j] * 0.25 *48 #Simultaneous conversion from density to quantity and summation \n            \n    return no2amount","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def calc_own_ef(df, no2amount):\n    \"\"\"\n    Calculate own emission factor of each power plant.\n    \n    Parameters\n    ----------\n    df : pandas.DataFrame\n        This dataframe must have gwh_rate_group and kmean_group.\n    no2amount : dictionary\n        value: group number, value: amount of substance of NO2\n\n    Returns\n    -------\n    df : pandas.DataFrame\n        Information of powerplant added estimated merginal_EF and emission for one year.\n    \"\"\"\n    own_efs = []\n    own_emission = []\n    \n    for i in range(len(df)):\n        pplant_cap_rate_group = df.iloc[[i]].loc[:,[\"cap_rate_group\"]].iat[0,0]\n        pplant_est_gen_gwh = df.iloc[[i]].loc[:,[\"estimated_generation_gwh\"]].iat[0,0]\n        pplant_kmean_group = df.iloc[[i]].loc[:,[\"kmean_group\"]].iat[0,0]\n        \n        #I calculate emissions factor of each power plant here!\n        #I assumed that one Sentinel-5P data coressponding to daily data.\n        own_efs.append(no2amount[pplant_kmean_group] * pplant_cap_rate_group / pplant_est_gen_gwh) \n        own_emission.append(no2amount[pplant_kmean_group] * pplant_cap_rate_group) \n        \n    df[\"own_EF\"] = own_efs\n    df[\"own_emission\"] = own_emission\n    \n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"no2amount_dict_group = calc_no2amount_each_group(image, group_pred)\n\nprint(\"Coreration of 'group: Total NO2 amount of each group(g)' are following:\")\nno2amount_dict_group #Total NO2 amount of each group","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"power_plants_df = calc_own_ef(power_plants_df, no2amount_dict_group)\npower_plants_df.loc[:,[\"name\", \"primary_fuel\", 'estimated_generation_gwh','capacity_mw' , \"kmean_group\", \"own_emission\", \"own_EF\"]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, we successed to get emissions factor of each power plant (own_EF)!\n\nNote that unit of own_EF is [g/GWh] amd own_emission is [g].\n\nAnd if the group and the fuel were the same, they would have the same value. This may be because capacity is used to calculate the estimated_generation_gwh\t. We think that it will be solved if we use the power generation results of last year instead of capacity."},{"metadata":{},"cell_type":"markdown","source":"## Find merginal emissions factor\n\nWe got emission factors of each power plants. Using these values, I'll calculate merginal emission factor.\n\nTo calculate merginal emission factor, we have to know following infomation, but we don't have them. \n\n- Power generation cost of each power plants\n\n- Which regions receive power from which power plants \n\nSo I only calculate merginal emission factor of group0 area assuming that the region of Group 0 depends on the power plants in the region."},{"metadata":{"trusted":false},"cell_type":"code","source":"power_plants_group0_df = power_plants_df[power_plants_df[\"kmean_group\"]==0]\npower_plants_group0_df.loc[:,[\"name\", \"primary_fuel\", \"kmean_group\", \"own_emission\", \"own_EF\"]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"total_capacity_mw_group0 = power_plants_group0_df['capacity_mw'].sum()\nprint('Total Installed Capacity: '+'{:.2f}'.format(total_capacity_mw_group0) + ' MW')\ncapacity = (power_plants_group0_df.groupby(['primary_fuel'])['capacity_mw'].sum()).to_frame()\ncapacity = capacity.sort_values('capacity_mw',ascending=False)\ncapacity['percentage_of_total'] = (capacity['capacity_mw']/total_capacity_mw_group0)*100\nax = capacity.sort_values(by='percentage_of_total', ascending=True)['percentage_of_total'].plot(kind='bar',color=['lightblue', 'green', 'orange', 'black','lightgray','darkblue'])\nax.set(ylabel='percentage')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"capacity","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"~9% of capacity, marginal emission factor is 0.\n\n9% ~ 41.5 of capacity, marginal emission factor is 0.071325 g/GWh.\n\n41.5% ~  of capacity, marginal emission factor is 0.089964 g/GWh."},{"metadata":{},"cell_type":"markdown","source":"# What is great of this method?"},{"metadata":{},"cell_type":"markdown","source":"## Pros and cons\n\n**<u>Pros</u>**\n\n- We can calculate emissions factor very simply. But the value seems reasonable. (method 1 and 2)\n- By inputting together weather and geographical data, it is possible to attribute which NOx is attributable to which power plant.(method 2)\n- Other useful data can be added to clustering simply by standardizing and adding.(method2)\n\n**<u>Cons</u>**\n\n- We need to decide the number of clusters ourselves. but the result of clustering desn't  always reflect our intention.\n- When we divide the data into clusters, if there are no power plants in that cluster, NO2 in that area will not be reflected in any power plants\n  (However, this disadvantage can be covered by distributing NO2 in these areas to all power plants.)\n- Data observed by remote sensing includes NOx from sources other than power plants such as viecle, bussiness site and so on. With this method, extra data is mixed into the emissions factor."},{"metadata":{},"cell_type":"markdown","source":"## Why this approach improve emissions factor?\n\n- We can calculate the NO2 emission factor directly from the actual amount present in the atmosphere. In \"Why pay attention to NO2?\", I explained that NOx reacts with chemicals in the atmosphere and changes into various forms. To know the NO2 density, it was necessary to separately investigate the concentration of ozone and the like and estimate the concentration that has reached equilibrium as a result of the chemical change.\n\n- To estimate which NO2 comes from which power plant, we can use weather, geographical data and so on. This time, I used only weather data, but if we find more efficient data, we can include them by just inputting to clustering model togather. "},{"metadata":{},"cell_type":"markdown","source":"## Reference\n\n[1] - [2] are documents by japanese public institutions. (\nSo sorry in Japanese.)\n\n[1]https://ghg-santeikohyo.env.go.jp/files/calc/cm_ec_R01/full.pdf\n\n[2] https://www.pref.chiba.lg.jp/taiki/shingikai/kankyou-taiki/documents/20110601shiryou04.pdf\n\n[3] - [5] are my kernel published in advance.\n\n[3]https://www.kaggle.com/nayuts/can-we-attribute-emissions-to-power-plants\n\n[4]https://www.kaggle.com/nayuts/focus-on-specific-power-plants\n\n[5]https://www.kaggle.com/nayuts/exploration-of-tiff-bands\n\n[6]https://www.kaggle.com/nayuts/calculate-ef-in-smaller-time-slices"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.6"}},"nbformat":4,"nbformat_minor":4}