{"cells":[{"metadata":{},"cell_type":"markdown","source":"## <u>About this kernel</u>\n\nTo calculate marginal emissions factor, I want to attribute emissions to each power plants.\n\n※For marginal emissions factor, refer https://www.tmrow.com/blog/marginal-emissions-what-they-are-and-when-to-use-them\n\n\nFor calculate reasonable emissions of each power plants, I try to apply k-means to satellite data as one of the typical clustering methods.\n\nIf we devide into small area, we can detect its area(m^2) and gas emission by following fomula:\n\n## Emissions(mol) = NO2_column_number_density(mol/m^2) * area(m^2)\n\nWhen we get mol of emited gas, we can calculate its weight (such as gram).\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"from decimal import Decimal, ROUND_HALF_UP, ROUND_HALF_EVEN\nimport folium\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":"## Snippts"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def overlay_image_on_puerto_rico_df(df, img, zoom):\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)\n    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\n    folium.raster_layers.ImageOverlay(\n        image=img,\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":true},"cell_type":"code","source":"def 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":{},"cell_type":"markdown","source":"## Data overview"},{"metadata":{"trusted":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() \npower_plants_df = power_plants.sort_values('capacity_mw',ascending=False).reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"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\ncapacity.sort_values(by='percentage_of_total', ascending=True)['percentage_of_total'].plot(kind='bar',color=['lightblue', 'green', 'orange', 'black','lightgray','darkblue'])\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Power plants and NO2_column_number_density of 20180701-20180707 in Puerto Rico**"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"image = 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#https://www.kaggle.com/paultimothymooney/explore-image-metadata-s5p-gfs-gldas\n#band1: NO2_column_number_density","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Especially, Coal, Gas and Oil power plants emit NO2 gas. So it may be nice idea that whole data divide into some areas. \n\nIdeally I want to devide whole data into areas as many as power plants in Puerto Rico, but it is difficult. One of the difficulty is that there are some power plants very near each other. \n\n\nIn small divided areas, it would be acceptable to divide the emissions by the amount of electricity generated."},{"metadata":{},"cell_type":"markdown","source":"# Devide NO2 map by k-means\n\nI try to generate devided areas by longitude, latitude and average NO2 density in a week."},{"metadata":{"trusted":true},"cell_type":"code","source":"lon = []\nlat = []\nNO2 = []\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        \nNO2 = np.array(NO2)\nlon = np.array(lon)\nlat = np.array(lat)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"results = pd.DataFrame(columns=['NO2', 'lat', 'lon'])\nresults = pd.DataFrame({'NO2': NO2/max(NO2),\n                    'lat': lat/max(lat),\n                    'lon': lon/max(lon)})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(results[\"NO2\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred = KMeans(n_clusters=11).fit_predict(results)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure()\nsns.heatmap(pred.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_df, pred.reshape((148, 475)), 8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Consideration\n\nConsidering the area made, there are some characteristics like below:\n\n・The devided areas seems be stronglly  constrained by coordinates. \n\n・ We have to decide the number of cluster manually.\n\nIt is nice that model more reasonablly devide while data into area if we give the data and coodinates of power plants.\n\nIf we input terrain and weather data to k-means, I think we may get more apposite attribute."},{"metadata":{},"cell_type":"markdown","source":"# Classification including wind data."},{"metadata":{},"cell_type":"markdown","source":"### Can we devide only by NO2 density?\n\nIn the first place, NO2 gas is distributes as we image? For example, if NO2 gas is distributed in a circle or a band from the power plant, it seems good to simply cluster at the density of NO2. But if not so, we have to use other data they seems efffect NO2 gas distribution.\n\nFitst, simply we can use k-means to convert no2 density to more monotonous classes. "},{"metadata":{"trusted":true},"cell_type":"code","source":"monotonous = KMeans(n_clusters=3).fit_predict(image[:,:,0].reshape(-1, 1))\n#pred = KMeans(n_clusters=2).fit_predict(results[\"NO2\"])\nplt.figure()\nsns.heatmap(monotonous.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Note that the color intensity and the high NO2 concentration do not always match.\noverlay_image_on_puerto_rico_df(power_plants_df, monotonous.reshape((148, 475)), 8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\nAs you can see from this figure, simply clustering seems to be difficult because the NO2 gas concentration is not distribute in a circle or a band from the power plant.　😅"},{"metadata":{},"cell_type":"markdown","source":"## Including Wind Data\n\nI tried to make wind feature, and make cluster by k-means.\n\nNO2 gases flowed by wind, and accumulate place where no wind."},{"metadata":{"trusted":true},"cell_type":"code","source":"import os \nimport glob","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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":"Let's retrieve the data from 7/1 to 7/7. And I put data together each day."},{"metadata":{"trusted":true},"cell_type":"code","source":"gldas_files_par_day = []\nfor i in range(0,len(gldas_files[6:54]),8):\n    #print(gldas_files[i:i+8])\n    gldas_files_par_day.append(gldas_files[i:i+8])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gfs_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":"### future engineering\n\nI tried flooowing 2 way.\n\nWe can use following three wind propaties: wind_u, wind_v and wind speed. They are given smaller time scale than NO2 gas.\n\n1. just calculate average of wind_u, wind_v and wind speed.\n\n2. make feature using wind_u, wind_v and wind speed. These are calculated by multipling wind speed to wind_u and wind_v."},{"metadata":{},"cell_type":"markdown","source":"### <u>Way1</u>"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_reglession_u = []\nimage_reglession_v = []\nimage_reglession_speed = []\n\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    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        \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        image_reglession_u.append(np.nanmean(np.array(array_wind_u), axis=0))\n        image_reglession_v.append(np.nanmean(np.array(array_wind_v), axis=0))\n        image_reglession_speed.append(np.nanmean(np.array(array_wind_speed), axis=0))\n       \nimage_reglession_u = np.nanmean(np.array(image_reglession_u), axis=0)\nimage_reglession_v = np.nanmean(np.array(image_reglession_v), axis=0)\nimage_reglession_speed = np.nanmean(np.array(image_reglession_speed), axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_u.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_v.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_speed.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = tiff.imread('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif')\nlon = []\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(image_reglession_u.reshape((148, 475))[i,j])\n        wind_v.append(image_reglession_v.reshape((148, 475))[i,j])\n        wind_speed.append(image_reglession_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        \nresults_wind = pd.DataFrame(columns=['NO2', 'lat', 'lon', 'wind_u', 'wind_v', 'wind_speed'])\nresults_wind = 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":{"trusted":true},"cell_type":"code","source":"pred_wind1 = KMeans(n_clusters=11).fit_predict(results_wind)\nplt.figure()\nsns.heatmap(pred_wind1.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_df, pred_wind1.reshape((148, 475)), 8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### <u>Way2</u>"},{"metadata":{"trusted":true},"cell_type":"code","source":"image_reglession_u = []\nimage_reglession_v = []\n\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    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        \n        gldas_image = (gldas_image1 + gldas_image2)/2\n        wind_u = gfs_image_u * gldas_image\n        wind_v = gfs_image_v * gldas_image\n        \n        array_wind_u.append(wind_u)\n        array_wind_v.append(wind_v)\n        \n        image_reglession_u.append(np.nanmean(np.array(array_wind_u), axis=0))\n        image_reglession_v.append(np.nanmean(np.array(array_wind_v), axis=0))\n       \nimage_reglession_u = np.nanmean(np.array(image_reglession_u), axis=0)\nimage_reglession_v = np.nanmean(np.array(image_reglession_v), axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_u.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(image_reglession_v.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = tiff.imread('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif')\nlon = []\nlat = []\nNO2 = []\nwind_u = []\nwind_v = []\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(image_reglession_u.reshape((148, 475))[i,j])\n        wind_v.append(image_reglession_v.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)\n        \nresults_wind = pd.DataFrame(columns=['NO2', 'lat', 'lon', 'wind_u', 'wind_v'])\nresults_wind = 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\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_wind2 = KMeans(n_clusters=11).fit_predict(results_wind)\nplt.figure()\nsns.heatmap(pred_wind2.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"overlay_image_on_puerto_rico_df(power_plants_df, pred_wind2.reshape((148, 475)), 8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Conparison"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure()\nsns.heatmap(pred.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure()\nsns.heatmap(pred_wind1.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure()\nsns.heatmap(pred_wind2.reshape((148, 475)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}