{"cells":[{"metadata":{},"cell_type":"markdown","source":"I want to calculate EF of each power plant, finally.\n\nSo, the procedure is\n    1. NO2 image data (tropo, strato, ... etc)\n    2. Weather(gfs) image data at band (temp, wind, ...etc)\n    3. To check NO2 density and wether day by day\n    4. To estimate to devide NO2 density caused by power plant and others\n        (Compare the NO2 density near plant and far.)\n    5. calculate EF of power plant"},{"metadata":{},"cell_type":"markdown","source":"Import liblary and csv file"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport os\nimport rasterio as rio\nimport folium\nimport tifffile as tiff \n\nimport seaborn as sns\n\nimport datetime as dt\nfrom datetime import datetime \n\nfrom decimal import Decimal, ROUND_HALF_UP, ROUND_HALF_EVEN\nimport math\n\nimport folium\nimport rasterio as rio\n\nimport warnings\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"global_power_plants = pd.read_csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Create NO2 Image data at band 1 to 7"},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_im(file,dir_path):\n    return rio.open(dir_path+file)\n\ndef mean_yearly_image(band):\n    dir_s5p_no2 = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/'\n    all_data = [read_im(file,dir_s5p_no2).read(band) for file in os.listdir(dir_s5p_no2)]\n    return np.array(all_data)\n\nall_im_no2 = mean_yearly_image(1)\nall_im_tropono2 = mean_yearly_image(2)\nall_im_stratono2 = mean_yearly_image(3)\nall_im_slantno2 = mean_yearly_image(4)\n\nall_im_troppresure = mean_yearly_image(5)\nall_im_absindex= mean_yearly_image(6)\nall_im_cloudfraction= mean_yearly_image(7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Upside Down NO2 image file\n##because fgs image and NO2 image is oposition latitude\n\nall_im_no2=np.flip(all_im_no2,axis=1)\nall_im_tropono2=np.flip(all_im_tropono2,axis=1)\nall_im_stratono2=np.flip(all_im_stratono2,axis=1)\nall_im_slantno2 = np.flip(all_im_slantno2,axis=1)\nall_im_troppresure =np.flip(all_im_troppresure,axis=1)\nall_im_absindex= np.flip(all_im_absindex,axis=1)\nall_im_cloudfraction= np.flip(all_im_cloudfraction,axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#NO2 density on Puerto Rico area\nsns.heatmap(np.nanmean(all_im_no2,axis=0),square=True,cmap='inferno',cbar_kws={'shrink':0.7})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.heatmap(np.nanmean(all_im_cloudfraction,axis=0),square=True,cmap='inferno',cbar_kws={'shrink':0.7})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Create gfs image data at band 1 ,3, 4, 5"},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_im(file,dir_path):\n    return rio.open(dir_path+file)\n\ndef mean_yearly_image(band):\n    dir_gfs = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/'\n    all_data = [read_im(file,dir_gfs).read(band) for file in os.listdir(dir_gfs)]\n    return np.array(all_data)\n\nall_im_temp = mean_yearly_image(1)\nall_im_humidity = mean_yearly_image(3)\nall_im_u_wind = mean_yearly_image(4)\nall_im_v_wind = mean_yearly_image(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Create Global Power Plant data\n#new column \"latitude\", \"longitude\"\n\nglobal_power_plants['latitude']=0\nglobal_power_plants['longitude']=0\nfor i in range(0, len(global_power_plants)):\n    global_power_plants['latitude'][i]=global_power_plants['.geo'][i][50:66]\n    global_power_plants['longitude'][i]=global_power_plants['.geo'][i][31:48]\n    \nglobal_power_plants['latitude'] = global_power_plants['latitude'].astype(float)\na = np.array(global_power_plants['latitude'].values.tolist()) # 18 instead of 8\nglobal_power_plants['latitude'] = np.where(a < 10, a+10, a).tolist() \n\n#Delete useless columns\n\npower_plants=global_power_plants.drop('generation_gwh_2013',axis=1)\npower_plants=power_plants.drop('generation_gwh_2014',axis=1)\npower_plants=power_plants.drop('generation_gwh_2015',axis=1)\npower_plants=power_plants.drop('generation_gwh_2016',axis=1)\npower_plants=power_plants.drop('generation_gwh_2017',axis=1)\npower_plants=power_plants.drop('other_fuel1',axis=1)\npower_plants=power_plants.drop('other_fuel2',axis=1)\npower_plants=power_plants.drop('other_fuel3',axis=1)\npower_plants=power_plants.drop('year_of_capacity_data',axis=1)\npower_plants=power_plants.drop('system:index',axis=1)\npower_plants=power_plants.drop('country',axis=1)\npower_plants=power_plants.drop('country_long',axis=1)\npower_plants=power_plants.drop('geolocation_source',axis=1)\npower_plants=power_plants.drop('gppd_idnr',axis=1)\npower_plants=power_plants.drop('url',axis=1)\npower_plants=power_plants.drop('wepp_id',axis=1)\npower_plants=power_plants.drop('.geo',axis=1)\n\npower_plants.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Transform longitude and latitude to pixel of tif file data array."},{"metadata":{"trusted":true},"cell_type":"code","source":"import rasterio as rio\n\ngfs_file = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/gfs_2018070206.tif'\n\ndef preview_meta_data(file_name):\n    with rio.open(file_name) as img_filename:\n        print('Bounding Box:',img_filename.bounds)\n        print('Coordinates of Top Left Corner: ',img_filename.transform * (0,0))\n        print('Coordinates of Bottom Right Corner: ',img_filename.transform * (img_filename.width,img_filename.height))\n        print(img_filename.index(-65.19081511310455, 18.564903861343627))\n        print(img_filename.index(-67.32354977311168, 18.564903861343627))\n        \npreview_meta_data(gfs_file)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#1pixel is equal to ...\nleft=-67.32354977311168\nright=-65.19005097332781\ntop=18.56520446703891\nbottom=17.900451156790464\n\npixel_length=(right-left)/474*225216.070 #[m]　Distance is calculated by https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2stf.html\npixel_width=(top-bottom)/147*73576.937 #[m]\n\npixel_area=pixel_length*pixel_width #[m^2]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_name = gfs_file\n\nwith rio.open(file_name) as img_filename:\n    power_plants['gfs_pixel_x']=img_filename.index(power_plants.longitude.astype(float), power_plants.latitude.astype(float))[0]\n    power_plants['gfs_pixel_y']=img_filename.index(power_plants.longitude.astype(float), power_plants.latitude.astype(float))[1]\n    \npower_plants","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To simplify, focus on Vieques EPP as example"},{"metadata":{"trusted":true},"cell_type":"code","source":"#Check vieques power plant data, again\nvieques_data=power_plants.loc[19:19]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Transform NO2 image data and gfs image data to dataframe day by day\n\nimport os\n\nno2_path='../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/'\nno2_list=os.listdir('../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2')\nno2_list.sort()\n\ngfs_path='../input/ds4g-environmental-insights-explorer/eie_data/gfs/'\ngfs_list=os.listdir('../input/ds4g-environmental-insights-explorer/eie_data/gfs/')\ngfs_list.sort()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get date and time from image file name\nfrom datetime import datetime, timedelta\ndatatime_no2=[datetime.strptime(i[:16],'s5p_no2_%Y%m%d') for i in no2_list]\ndatatime_gfs=[datetime.strptime(i,'gfs_%Y%m%d%H.tif') for i in gfs_list]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Make daily data of Vieques Power Plant at fgs_pixel_x and y\n#To compare the NO2 density far from power plant, make data at fgs pixel_x+10,y+10\n#(because maily oposition direction of wind)\n\nvieques_daily_no2=pd.DataFrame(datatime_no2)\nvieques_daily_gfs=pd.DataFrame(datatime_gfs)\n\n#Mean of NO2 density around Vieques Power Plant\nvieques_daily_no2['tropono2']=all_im_tropono2[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_no2['no2']=all_im_no2[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\n\nvieques_daily_no2['no2_max']=all_im_no2[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_no2['no2_min']=all_im_no2[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\n\nfor k in range(len(all_im_no2)):\n    for i in range(-5,5):\n        for j in range(-5,5):\n            vieques_daily_no2.no2_max[k]=max(vieques_daily_no2.no2_max[k],all_im_no2[k,vieques_data.gfs_pixel_x+i,vieques_data.gfs_pixel_y+j])\n            vieques_daily_no2.no2_min[k]=min(vieques_daily_no2.no2_min[k],all_im_no2[k,vieques_data.gfs_pixel_x+i,vieques_data.gfs_pixel_y+j])\n\nvieques_daily_no2['stratono2']=all_im_stratono2[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_no2['slantno2']=all_im_slantno2[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_no2['troppresure']=all_im_troppresure[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_no2['absindex']=all_im_absindex[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_no2['cloudfraction']=all_im_cloudfraction[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\n\n#weather data at Vieques Power Plant\nvieques_daily_gfs['temp']=all_im_temp[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_gfs['humidity']=all_im_humidity[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_gfs['u_wind']=all_im_u_wind[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\nvieques_daily_gfs['v_wind']=all_im_v_wind[:,vieques_data.gfs_pixel_x,vieques_data.gfs_pixel_y]\n\n#Change column name to set_index\nvieques_daily_gfs=vieques_daily_gfs.rename(columns={0:'gfsdate'})\nvieques_daily_no2=vieques_daily_no2.rename(columns={0:'no2date'})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_daily_no2.set_index('no2date',inplace=True)\nvieques_daily_gfs.set_index('gfsdate',inplace=True)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_daily_no2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_daily_gfs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#There are four wether data each day(0:00,6:00,12:00,18:00), so select 12:00 data\nvieques_daily_gfs=vieques_daily_gfs[vieques_daily_gfs.index.hour == 12]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_daily_gfs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#To merge NO2 density data and weather data\nvieques_daily_gfs['day']=vieques_daily_gfs.index\nvieques_daily_no2['day']=vieques_daily_no2.index\n\nvieques_daily_gfs['day']=vieques_daily_gfs.day.dt.date\n\nvieques_daily_no2.day=vieques_daily_no2.day.astype('str')\nvieques_daily_gfs.day=vieques_daily_gfs.day.astype('str')\n\nvieques_daily=pd.merge(vieques_daily_no2,vieques_daily_gfs,on=['day'],how='inner')\nvieques_daily=vieques_daily[['day','no2','no2_max','no2_min','tropono2','stratono2', 'slantno2', 'troppresure',\n                             'absindex','cloudfraction','temp','humidity','u_wind','v_wind']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Make wind velocity data from u and v\nvieques_daily['wind']=np.sqrt(vieques_daily.u_wind**2+vieques_daily.v_wind**2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"### Final dataframe for Vieques Power plant\nvieques_daily","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 10))\nsns.distplot(vieques_daily.no2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see the correlation. To hope that there is correlation between NO2 and wind."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(vieques_daily)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I can't find clear correlation...."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 10))\nsns.jointplot('no2','temp',data=vieques_daily,xlim=(-0.0003,0.0003))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 10))\nsns.jointplot('no2','wind',data=vieques_daily,xlim=(-0.00003,0.0001))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To see the trend of NO2 desity, wind, "},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfig = plt.figure(figsize=(18, 6))\nplt.plot(vieques_daily.no2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfig = plt.figure(figsize=(18, 6))\nplt.plot(vieques_daily.wind)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfig = plt.figure(figsize=(18, 6))\nplt.plot(vieques_daily.absindex)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_daily.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To devide NO2 emitted from plant and other, calculate the difference NO2 density near and far from power plant."},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_daily['dif_no2']=vieques_daily.no2_max-vieques_daily.no2_min\nvieques_daily","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfig = plt.figure(figsize=(20, 6))\nplt.plot(vieques_daily.day,vieques_daily.dif_no2)\nplt.xlim(210,220)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#To sumarize the daily data to the monthly data, 'day' object change to datetime\nvieques_daily['day_pd']=pd.to_datetime(vieques_daily.day,format='%Y-%m-%d')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_daily.set_index('day_pd',inplace=True)\nvieques_daily","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_month=vieques_daily.resample(rule=\"M\").sum()\nvieques_month","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(vieques_month)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"NO2 has correlation with 'troppresure','temp' and 'humidity'.\nSo, it may be possible to estimate base density of NO2. (\"base\" means \"without power plant\".)\n\nBut now I can't do that, because of my skill..."},{"metadata":{},"cell_type":"markdown","source":"**To calculate emission** by difference of NO2 and 1 pixel area. Unit is ton, and NO2 is 46.0055g/mol.**"},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_month['emission']=vieques_month.dif_no2*pixel_area*46.0055/1000000\nvieques_month","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next , I will estimate generation yearly, and each month."},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants.loc[19:19]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"EF_year = sum(vieques_month.emission)/power_plants.loc[19:19].estimated_generation_gwh/24/365\nprint(\"EF_year = {}\" .format(EF_year[19])) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"vieques_month['EF_month']=EF_year[19]*vieques_month.emission/sum(vieques_month.emission)\nvieques_month","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}