{"cells":[{"metadata":{},"cell_type":"markdown","source":"# PROBLEM STATEMENT"},{"metadata":{},"cell_type":"markdown","source":"Current emissions factors methodologies are based on time-consuming data collection and may include errors derived from a lack of access to granular datasets, inability to refresh data on a frequent basis, overly general modeling assumptions, and inaccurate reporting of emissions sources like fuel consumption. This begs the question: What if there was a different way to calculate or measure emissions factors? We’re challenging the Kaggle community to see if it’s possible to use remote sensing techniques to better model emissions factors. You will develop a methodology to calculate an average historical emissions factor for electricity generation in a sub-national region.\n\nWe’ve provided an initial list of datasets covering the geographic boundary of Puerto Rico to serve as the foundation for this analysis. As an island, there are fewer confounding factors from nearby areas. Puerto Rico also offers a unique fuel mix and distinctive energy system layout that should make it easier to isolate pollution attributable to power generation in the remote sensing data.\n\nParticipants will be tasked with developing a methodology to calculate an average annual historical emissions factor for the sub-national region. Participants will also be asked to provide an explanation of the conditions that would result in a higher/lower emissions factor, as well as a recommendation for how the methodology could be applied to calculate the emissions factor of electricity for another geospatial area using similar techniques. Bonus points will be awarded for smaller time slices of the average historical emissions factors, such as one per month for the 12-month period, and additional bonus points will be awarded for participants that develop methodologies for calculating marginal emissions factors for the sub-national region."},{"metadata":{},"cell_type":"markdown","source":"# EXPLORATORY DATA ANALYSIS"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport rasterio as rio\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport os\nimport folium\nimport tifffile as tiff\nfrom datetime import datetime, timedelta\nfrom branca.element import Template, MacroElement\nfrom skimage.transform import resize\nimport branca\nfrom tqdm import tqdm\nfrom scipy.stats import pearsonr\nimport geopy.distance\nfrom sklearn.preprocessing import normalize\nimport cv2\nimport warnings\nfrom affine import Affine\nfrom pyproj import Proj, transform","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"warnings.simplefilter('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"According to the [American Lung Association](https://www.lung.org/clean-air/outdoors/what-makes-air-unhealthy/nitrogen-dioxide), cars, trucks, and buses are the largest sources of NOx emissions, followed by power plants, diesel-powered heavy construction equipment and other movable engines, and industrial boilers. Man-made sources in the U.S. emitted 14 million metric tons of nitrogen oxides, mainly from burning fuels, in 2011. Emissions of nitrogen dioxide will decline as cleanup of many of these sources continue in future years."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"s5p_no2_path = '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/'\ngfs_path = '../input/ds4g-environmental-insights-explorer/eie_data/gfs/'\n\nelectricity = pd.read_csv('../input/ds4g-puerto-rico-electricity-consumption-by-month/EIA_puerto_rico_electric_consumption.csv')\nmonth_av_elect = electricity.iloc[30:42]['all_sectors'].values\n#sorting from Jan to Dec\nmonth_av_elec = np.append(month_av_elect[6:],month_av_elect[:6])\n\ngppd = pd.read_csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\ngppd = gppd.drop(['system:index', 'capacity_mw', 'country',\n 'country_long', 'estimated_generation_gwh', 'generation_gwh_2013',\n 'generation_gwh_2014', 'generation_gwh_2015', 'generation_gwh_2016',\n 'generation_gwh_2017', 'geolocation_source', 'gppd_idnr',\n 'source', 'url', 'wepp_id', 'year_of_capacity_data'], axis=1)\ngppd['long'] = gppd['.geo'].apply(lambda x: float(x[31:48]))\ngppd['lat'] = gppd['.geo'].apply(lambda x: float(x[50:66]))\ngppd['lat'] = np.where(gppd['lat']<10, gppd['lat']+10, gppd['lat'])\ngppd['Ix']=gppd.index\ngppd['coords'] = list(zip(gppd.lat, gppd.long))\n\npop_map_file = '/kaggle/input/population/imageToDriveExample.tif'\npop_map = rio.open(pop_map_file).read(1)\n\neia = pd.read_csv('../input/eia923/PR_EIA_923.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"s5p_bounds = rio.open(s5p_no2_path+os.listdir(s5p_no2_path)[0]).bounds\ngfs_bounds = rio.open(gfs_path+os.listdir(gfs_path)[0]).bounds\ns5p_top, s5p_bot, s5p_lef, s5p_rit = s5p_bounds[1], s5p_bounds[3], s5p_bounds[0], s5p_bounds[2] #will be flipped while reading\ngfs_top, gfs_bot, gfs_lef, gfs_rit = gfs_bounds[3], gfs_bounds[1], gfs_bounds[0], gfs_bounds[2]\n\npop_bounds = rio.open(pop_map_file).bounds\npop_top, pop_bot, pop_lef, pop_rit = pop_bounds[3], pop_bounds[1], pop_bounds[0], pop_bounds[2]\n\n#get the x,y location on map image for each longitude, latitude\nno2_latp = (gppd['lat']-s5p_top)/(s5p_bot-s5p_top)*100\nno2_longp = (gppd['long']-s5p_lef)/(s5p_rit-s5p_lef)*100\nno2_col = np.round(475*no2_longp/100).astype(int)\nno2_row = np.round(148*no2_latp/100).astype(int)\n\npop_latp = (gppd['lat']-pop_top)/(pop_bot-pop_top)*100\npop_longp = (gppd['long']-pop_lef)/(pop_rit-pop_lef)*100\npop_col = np.round(256*pop_longp/100).astype(int)\npop_row = np.round(80*pop_latp/100).astype(int)\n\nwith rio.open(pop_map_file) as r:\n T0 = r.transform # upper-left pixel corner affine transform\n p1 = Proj(r.crs)\n A = r.read() # pixel values\n\ncols, rows = np.meshgrid(np.arange(A.shape[2]), np.arange(A.shape[1]))\n\nT1 = T0 * Affine.translation(0.5, 0.5)\n# Function to convert pixel row/column index (from 0) to easting/northing at centre\nrc2en = lambda r, c: (c, r) * T1\n\n# All eastings and northings\neastings, northings = np.vectorize(rc2en, otypes=[np.float, np.float])(rows, cols)\n\n# Project all longitudes, latitudes\np2 = Proj(proj='latlong',datum='WGS84')\nlongs, lats = transform(p1, p2, eastings, northings)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#color dictionary for map and plots\nfuels = ['Coal','Gas','Hydro','Oil','Solar','Wind']\ncolors = ['black','gray','green','darkred','orange','blue']\nfudic = {}\nfor k,v in zip(fuels, colors):\n fudic[k]=v","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#no2 columns grouped by day of the week\nweek_no2 = {}\nweek_no2_trop = {}\nweek_no2_strat = {}\nfor i in range(7):\n week_no2[i]=[]\n week_no2_trop[i] = []\n week_no2_strat[i] = []\nfor file in tqdm(os.listdir(s5p_no2_path)):\n ao = datetime.strptime(file[8:16], '%Y%m%d')\n filr = rio.open(s5p_no2_path+file)\n cf = np.flipud(filr.read(7))\n no = np.flipud(filr.read(1))\n tno = np.flipud(filr.read(2))\n sno = np.flipud(filr.read(3))\n #no = np.where(cf>0.3,np.nan,no)\n #tno = np.where(cf>0.3,np.nan,tno)\n #sno = np.where(cf>0.3,np.nan,sno)\n week_no2[ao.weekday()].append(no)\n week_no2_trop[ao.weekday()].append(no)\n week_no2_strat[ao.weekday()].append(no)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"all_yr_no2 = np.empty([0,148,475])\nfor i in range(7):\n all_yr_no2 = np.concatenate((all_yr_no2, np.array(week_no2[i])), axis=0)\nannual_no2_avg = np.nanmean(all_yr_no2, axis=0)\ncmap = (annual_no2_avg-annual_no2_avg.min())/(annual_no2_avg.max()-annual_no2_avg.min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gppd['annual_no2_av']=np.nan\ngppd['population_dens']=np.nan\nfor i in range(len(gppd)):\n plant=gppd.iloc[i]\n gppd.loc[i,'annual_no2_av']=annual_no2_avg[int(no2_row[i]),int(no2_col[i])]\n gppd.loc[i,'population_dens']=np.nanmean(pop_map[int(pop_row[i])-10:int(pop_row[i])+10,int(pop_col[i])-10:int(pop_col[i])+10])\n#gppd['population_dens']=gppd['population_dens'].fillna(0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# IMPACT OF POPULATION"},{"metadata":{},"cell_type":"markdown","source":"Let us take a look at the location of the power plants in relation to population centers, where vehicular exhaust is a major source of NOx."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gppd.groupby('primary_fuel')[['annual_no2_av','population_dens']].mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see above, the NO2 concentration depends on the type of power plant as well as the population density."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"popn =normalize(np.nan_to_num(pop_map))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"lat=(s5p_top+s5p_bot)/2; lon=(s5p_lef+s5p_rit)/2\nm = folium.Map([lat, lon], zoom_start=9, height=350, width = 1000)\nfolium.raster_layers.ImageOverlay(\n image=popn,\n bounds = [[pop_bot,pop_lef],[pop_top,pop_rit]],\n colormap=lambda x: (1, 0, 1, x),\n ).add_to(m)\nfor i in range(len(gppd)):\n folium.Marker([gppd['lat'].iloc[i],gppd['long'].iloc[i]],icon=folium.Icon(color=fudic[gppd['primary_fuel'].iloc[i]])).add_to(m)\ntitle_html = '''\n
Power Plants in Puerto Rico with Population Density
\n '''\nm.get_root().html.add_child(folium.Element(title_html))\ncolormap.add_to(m)\nm.get_root().add_child(macro)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"p = resize(pop_map, (148,475))\nfig, ax = plt.subplots(1,2, figsize=(18,12))\nax[0].axhline(y=0, color='k',linewidth=1)\nax[0].axhline(y=147, color='k',linewidth=1)\nax[0].axvline(x=0, color='k',linewidth=1)\nax[0].axvline(x=474, color='k',linewidth=1)\nsns.heatmap(p,ax=ax[0],square=True,cmap=sns.color_palette(\"Blues\"), cbar_kws={\"orientation\": \"horizontal\",'shrink':0.5})\nsns.heatmap(cmap, ax=ax[1], square=True, cmap='Reds',cbar_kws={\"orientation\": \"horizontal\",'shrink':0.5})\nsns.heatmap(p,ax=ax[1],square=True,cmap=sns.color_palette(\"Blues\"), cbar=False, alpha=0.2)\nfor fuel in fuels:\n i = gppd[gppd['primary_fuel']==fuel].index\n plt.scatter(no2_col[i], no2_row[i], color=fudic[fuel], label=fuel) \n#ax[0].axis('off')\n#ax[1].axis('off')\nax[0].set_title('Population Density')\nax[1].set_title('Relative NO2 Concentration with Population Density')\n#ax.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Densely populated areas around San Juan show the highest NOx concentration."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"weekend = np.concatenate((np.array(week_no2[5]),np.array(week_no2[6])), axis=0)\nweekday = np.concatenate((np.array(week_no2[0]),np.array(week_no2[1]),np.array(week_no2[2]),np.array(week_no2[3]),np.array(week_no2[4])), axis=0)\nweekend_trop = np.concatenate((np.array(week_no2_trop[5]),np.array(week_no2_trop[6])), axis=0)\nweekday_trop = np.concatenate((np.array(week_no2_trop[0]),np.array(week_no2_trop[1]),np.array(week_no2_trop[2]),np.array(week_no2_trop[3]),np.array(week_no2_trop[4])), axis=0)\nweekend_strat = np.concatenate((np.array(week_no2_strat[5]),np.array(week_no2_strat[6])), axis=0)\nweekday_strat = np.concatenate((np.array(week_no2_strat[0]),np.array(week_no2_strat[1]),np.array(week_no2_strat[2]),np.array(week_no2_strat[3]),np.array(week_no2_strat[4])), axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"weekendav = np.nanmean(weekend, axis=0)\nweekdayav = np.nanmean(weekday, axis=0)\nweekendav = np.where(np.isnan(p), np.nan,weekendav)\nweekdayav = np.where(np.isnan(p), np.nan,weekdayav)\n\nweekendav_trop = np.nanmean(weekend_trop, axis=0)\nweekdayav_trop = np.nanmean(weekday_trop, axis=0)\nweekendav_trop = np.where(np.isnan(p), np.nan,weekendav_trop)\nweekdayav_trop = np.where(np.isnan(p), np.nan,weekdayav_trop)\n\nweekendav_strat = np.nanmean(weekend_strat, axis=0)\nweekdayav_strat = np.nanmean(weekday_strat, axis=0)\nweekendav_strat = np.where(np.isnan(p), np.nan,weekendav_strat)\nweekdayav_strat = np.where(np.isnan(p), np.nan,weekdayav_strat)\n\nfig, ax = plt.subplots(1,2, figsize=(16,10))\nax = ax.ravel()\nsns.heatmap(weekdayav,ax=ax[0],square=True,cmap='inferno',cbar_kws={\"orientation\": \"horizontal\"}, vmax=6.518138194350083e-05, vmin=4.5764838692296064e-05)\nsns.heatmap(weekendav,ax=ax[1],square=True,cmap='inferno',cbar_kws={\"orientation\": \"horizontal\"}, vmax=6.518138194350083e-05, vmin=4.5764838692296064e-05)\n[x.axis('off') for x in ax]\n\nax[0].set_title('Weekday Total NO2')\nax[1].set_title('Weekend Total NO2')\n\nfig, ax = plt.subplots(1,2, figsize=(16,10))\nax = ax.ravel()\nsns.heatmap(weekdayav_trop,ax=ax[0],square=True,cmap='inferno',cbar_kws={\"orientation\": \"horizontal\"}, vmax=np.nanmax(weekdayav_trop), vmin=np.nanmin(weekendav_trop))\nsns.heatmap(weekendav_trop,ax=ax[1],square=True,cmap='inferno',cbar_kws={\"orientation\": \"horizontal\"}, vmax=np.nanmax(weekdayav_trop), vmin=np.nanmin(weekendav_trop))\n[x.axis('off') for x in ax]\n\nax[0].set_title('Weekday Tropospheric NO2')\nax[1].set_title('Weekend Tropospheric NO2')\n\nfig, ax = plt.subplots(1,2, figsize=(16,10))\nax = ax.ravel()\nsns.heatmap(weekdayav_strat,ax=ax[0],square=True,cmap='inferno',cbar_kws={\"orientation\": \"horizontal\"}, vmax=np.nanmax(weekdayav_strat), vmin=np.nanmax(weekendav_strat))\nsns.heatmap(weekendav_strat,ax=ax[1],square=True,cmap='inferno',cbar_kws={\"orientation\": \"horizontal\"}, vmax=np.nanmax(weekdayav_strat), vmin=np.nanmax(weekendav_strat))\n[x.axis('off') for x in ax]\n\nax[0].set_title('Weekday Stratospheric NO2')\nax[1].set_title('Weekend Stratospheric NO2')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The impact of motor vehicles can be seen in the difference in NOx profiles measured on weekdays vs weekends."},{"metadata":{},"cell_type":"markdown","source":"In a separate [notebook](https://www.kaggle.com/parselt/cloud-fraction/), we can see that the cloud fraction greatly influences the measurement of NO2, and needs to be accounted for in our analysis. It is also to be noted that as Puerto Rico is humid and cloudy all year round, excluding data points with high cloud fraction (>0.2-0.3) that has been done in other studies ([1](https://www.sciencedirect.com/science/article/abs/pii/S1352231015301291)) is not a possibility here. However, for areas with low cloud cover, these should be excluded."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"##nox dataframe of dates with corresponding files\nnodf = {}\nfor file in os.listdir(s5p_no2_path):\n ao = datetime.strptime(file[8:16], '%Y%m%d')\n nodf[ao] = file\n \nnodf = pd.DataFrame.from_dict(nodf, orient='index', columns=['file'])\nnodf['date'] = nodf.index\nnodf = nodf.sort_values(by='date')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"month_no2 = {}\nmonth_cf = {}\nfor mth in range(1,13):\n month_no2[mth]=[]\n month_cf[mth]=[]\nfor file in nodf.file:\n ao = datetime.strptime(file[8:16], '%Y%m%d')\n filr = rio.open(s5p_no2_path+file)\n no = np.flipud(filr.read(1))\n cf = np.flipud(filr.read(7))\n nocf = np.divide(no, 1-0.95*cf)\n nocf = np.nan_to_num(nocf, posinf=np.nan)\n month_no2[ao.month].append(nocf)\n month_cf[ao.month].append(cf)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"wind, angle, temp, humidity = {}, {}, {}, {}\nfor mth in range(1,13):\n wind[mth]=[]\n angle[mth]=[]\n temp[mth]=[]\n humidity[mth]=[]\nfor dt in tqdm(pd.to_datetime(nodf.date)):\n file = gfs_path + 'gfs_'+ dt.strftime('%Y%m%d') + '12.tif'\n gfs_im = rio.open(file)\n us = gfs_im.read(4)\n vs = gfs_im.read(5)\n wind[dt.month].append(np.sqrt(us**2+vs**2))\n angle[dt.month].append(np.arctan2(vs,us))\n temp[dt.month].append(gfs_im.read(1))\n humidity[dt.month].append(gfs_im.read(2))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Please note that months January through June are in 2019, and the rest are in 2018."},{"metadata":{"trusted":true},"cell_type":"code","source":"month_av_no2 ={}\nmonth_av_cf ={}\nmonth_av_wind = {}\nmonth_av_ang = {}\nmonth_av_temp ={}\nmonth_av_hum ={}\nfor i in range(1,13):\n month_av_no2[i] = np.nanmean(np.array(month_no2[i]), axis=0)\n month_av_cf[i] = np.nanmean(np.array(month_cf[i]), axis=0)\n month_av_wind[i] = np.nanmean(np.array(wind[i]), axis=0)\n month_av_ang[i] = np.nanmean(np.array(angle[i]), axis=0)\n month_av_temp[i] = np.nanmean(np.array(temp[i]), axis=0)\n month_av_hum[i] = np.nanmean(np.array(humidity[i]), axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"avno = []\navcf = []\navwind = []\navtemp = []\navhum = []\nfor i in range(1,13):\n avno.append(np.nanmean(month_av_no2[i]))\n avcf.append(np.nanmean(month_av_cf[i]))\n avwind.append(np.nanmean(month_av_wind[i]))\n avtemp.append(np.nanmean(month_av_temp[i]))\n avhum.append(np.nanmean(month_av_hum[i]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3, figsize=(20,5))\nax[0].plot(avno, label='NOx')\nax[0].twinx().plot(avhum, label='Humid', color='red')\nax[1].plot(avwind, label='Wind', color='green')\nax[2].plot(avtemp, label='Temp', color='yellow')\nax[2].twinx().plot(avcf, label='CF', color='orange')\nfig.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# ANNUAL EMISSIONS FACTOR"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,1)\nplt.plot(avno, label='Avg NOx')\nax2 = ax.twinx()\nplt.plot(month_av_elec, color='red', label='Total Electricity Gen')\nfig.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Emitted NOx closely correlates with electricity generation for the entire island of Puerto Rico."},{"metadata":{"trusted":true},"cell_type":"code","source":"month_ef = (avno/month_av_elec)*10**8 #multiply 10^8 for convenience\nmonth_ef.mean(), month_ef.std()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.kdeplot(month_ef)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Taking the average of the monthly emissions factors, we estimate a monthly emissions factor of 7.15 for Puerto Rico. This assumes that the proportion of NO2 generated from power plants to total NO2 generation in Puerto Rico is a constant over months. There maybe seasonal variations in traffic but in the absence of data, this is a reasonable assumption."},{"metadata":{},"cell_type":"markdown","source":"# INDIVIDUAL POWER PLANTS"},{"metadata":{},"cell_type":"markdown","source":"We now take a look at individual power plants to see if we can estimate the emissions factors at source level. Monthly electricity generation is available for individual plants in 2018 from [EIA](https://www.eia.gov/electricity/data/eia923/)."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"mths = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']\nno2cols = [mth+'_no2' for mth in mths]\ntempcols = [mth+'_temp' for mth in mths]\nhumidcols = [mth+'_hum' for mth in mths]\nwindcols = [mth+'_wind' for mth in mths]\nangcols = [mth+'_ang' for mth in mths]\nfor col in no2cols:\n gppd[col]=np.nan\nfor col in tempcols:\n gppd[col]=np.nan\nfor col in humidcols:\n gppd[col]=np.nan\nfor col in windcols:\n gppd[col]=np.nan\nfor col in angcols:\n gppd[col]=np.nan","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"for plant in range(len(gppd)):\n plantrow = gppd.iloc[plant]\n for i in range(12):\n gppd.loc[plant,mths[i]+'_no2']=month_av_no2[i+1][int(no2_row[plant]),int(no2_col[plant])]\n gppd.loc[plant,mths[i]+'_temp']=month_av_temp[i+1][int(no2_row[plant]),int(no2_col[plant])]\n gppd.loc[plant,mths[i]+'_hum']=month_av_hum[i+1][int(no2_row[plant]),int(no2_col[plant])]\n gppd.loc[plant,mths[i]+'_wind']=month_av_wind[i+1][int(no2_row[plant]),int(no2_col[plant])]\n gppd.loc[plant,mths[i]+'_ang']=month_av_ang[i+1][int(no2_row[plant]),int(no2_col[plant])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"nd=2\nfor plant in range(len(gppd)):\n plantrow = gppd.iloc[plant]\n for i in range(12):\n gppd.loc[plant,mths[i]+'_no2']=np.nanmean(month_av_no2[i+1][no2_row[i]-nd:no2_row[i]+nd,no2_col[i]-nd:no2_col[i]+nd])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#distance matrix to find the nearest plant and its primary fuel type for each power plant.\ndist = np.zeros((35,35))\n#dist[:] = np.nan\nfor i in range(len(gppd)):\n for j in range(i+1,len(gppd)):\n dist[i,j] = geopy.distance.distance(gppd['coords'].values[i],gppd['coords'].values[j]).km\ndist = dist + dist.T - np.diag(np.diag(dist))\nnear_dist = []\nnear_type = []\nnear_ix = []\nfor i in range(35):\n near_index = np.argsort(dist[i])[1]\n near_dist.append(dist[i][near_index])\n near_type.append(gppd['primary_fuel'].values[near_index])\n near_ix.append(near_index)\ngppd['near_dist']=near_dist\ngppd['near_type']=near_type\ngppd['near_index']=near_ix","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gppd[['name','near_dist','near_index','primary_fuel']].sort_values(by='near_dist')[:25]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some NOx generating plants are located in the same complex or very close to other power plants, so the plumes cannot be distinguished at the 5.5km x 3.5km resolution of the satellite. We can club these together for the sake of simplicity. The groups of relevance are -\n\n1. Palo Seco (17) and Palo Seco (3)\n2. Aguirre (1) and Aguirre (10)\n3. Costa Sur (2) and Costa Sur (12) with EcoElectrica (14)\n4. AES (28) and Jobos (15) "},{"metadata":{"trusted":true},"cell_type":"code","source":"gppd2 = pd.merge(gppd, eia, on='Ix', how='left')\nnetgencols = [x for x in gppd2.columns if 'Netgen' in x]\nfor ncol in netgencols:\n gppd2[ncol]=np.where(gppd2[ncol]=='.',0,gppd2[ncol])\n gppd2[ncol] = gppd2[ncol].str.replace(',', '').astype(float)\n gppd2[ncol]=np.where(gppd2[ncol].isna(),0,gppd2[ncol])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Only coal, oil and gas plants which generate NOx from burning fossil fuels are of interest to us."},{"metadata":{"trusted":true},"cell_type":"code","source":"gppdf = gppd2[gppd2['primary_fuel'].isin(['Coal','Oil','Gas'])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sap = ['Yabucoa', 'Daguao','Mayagüez', 'Vega Baja', 'Vieques EPP','Cambalache', 'San Juan CC']\nfig, ax = plt.subplots(3,3, figsize=(20,15))\n\nfor i in range(7):\n g = gppdf[gppdf['name']==sap[i]]\n r = int((i/3)%3)\n c = int(i%3)\n ax[r,c].plot(g[netgencols].sum()[6:], label='elec')\n ax[r,c].twinx().plot(g[no2cols].sum()[6:], color='red', label='no2')\n ax[r,c].set_title(sap[i])\nfig.legend()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gppdf1 = gppdf[gppdf['name']=='Aguirre']\ngppdf2 = gppdf[gppdf['name']=='Palo Seco']\ngppdf3 = gppdf[gppdf['name'].isin(['EcoEléctrica','Costa Sur'])]\ngppdf4 = gppdf[gppdf['name'].isin(['Jobos','A.E.S. Corp.'])]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,2, figsize=(20,5))\nax[0,0].plot(gppdf1[netgencols].sum()[6:], label='elec')\nax2=ax[0,0].twinx()\nax2.plot(gppdf1[no2cols].sum()[6:], color='red', label='no2')\n\nax[0,1].plot(gppdf2[netgencols].sum()[6:])\nax3=ax[0,1].twinx()\nax3.plot(gppdf2[no2cols].sum()[6:], color='red')\n\nax[1,0].plot(gppdf3[netgencols].sum()[6:])\nax4=ax[1,0].twinx()\nax4.plot(gppdf3[no2cols].sum()[6:], color='red')\n\nax[1,1].plot(gppdf4[netgencols].sum()[6:])\nax5=ax[1,1].twinx()\nax5.plot(gppdf4[no2cols].sum()[6:], color='red')\n\nfig.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plant level data is limited and may not be accurate. For more granularity in the emission factor, there is a need for verified plant level electricity generation data over a longer timeline. Atleast a few years of data is needed to account for seasonal variation. For plants in cloudless areas, the NO2 concentration would follow a Gaussian around the source. A 2-d wind adjusted Gaussian can be fitted to the measured NO2 data to estimate the source emission levels."}],"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}