{"cells":[{"metadata":{},"cell_type":"markdown","source":"**PORJECT OVERVIEW**\n\n***Develop a methodology to calculate an average historical emissions factor of electricity generated for a sub-national region, using remote sensing data and techniques.***"},{"metadata":{},"cell_type":"markdown","source":"**PROBLEM STATEMENT**\n\n***Measuring Emissions factors from Satellite Data***\n\n**To develop a methodology to calculate an average historical emissions factor for electricity generation in a sub-national region. **"},{"metadata":{},"cell_type":"markdown","source":"**DATA PROVIDED**\n\nInitial 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. Calculate the annual historical emission factor\n\nBonus points will be awarded for smaller time slices of the average historical emissions factors, such as one per month for the 12-month period\nAdditional 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":"* An emissions factor is a representative value that attempts to relate the quantity of a pollutant released to the atmosphere with an activity associated with the release of that pollutant. \n\n\n> These factors are usually expressed as the weight of pollutant divided by a unit weight, volume, distance, or duration of the activity emitting the pollutant (e.g., kilograms of particulate emitted per megagram of coal burned).\n\n> Such factors facilitate estimation of emissions from various sources of air pollution. In most cases, these factors are simply averages of all available data of acceptable quality, and are generally assumed to be representative of long-term averages for all facilities in the source category (i.e., a population average)."},{"metadata":{},"cell_type":"markdown","source":"**The general equation for emissions estimation is:**\n\n*E = A x EF x (1-ER/100)*\n\n**where:**\n\n    E = emissions;\n    \n    A = activity rate;\n    \n    EF = emission factor, and\n    \n    ER =overall emission reduction efficiency, %\n\n\n\n"},{"metadata":{},"cell_type":"markdown","source":"**PROBLEM APPROACH OR METHODOLOGY**\n\n    1. Analyse the temperature, wind speed and wind direction data to see if they correlate with the NO2 concentration \n    2. Perform principal component analysis to isolate pollution attributable to the power units\n    3  calculate the emission of NO2\n    4. calculate the emission factor\n    5. calculate the marginal emission factor\n    \n  "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport os\nimport rasterio as rio\nimport folium \nimport seaborn as sns\nfrom kaggle_secrets import UserSecretsClient\nfrom google.oauth2.credentials import Credentials\nfrom os import listdir\nimport ee\n\n# Import the Image function from the IPython.display module. \nfrom IPython.display import Image\nfrom folium import plugins \nfrom PIL import Image\n# Analysing datetime\nimport datetime as dt\nfrom datetime import datetime \n\nimport cartopy.crs as ccrs\nfrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\nimport matplotlib.ticker as mticker\n\n\n%matplotlib inline\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, t, filenames in os.walk('/kaggle/input'): \n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# Any results you write to the current directory are saved as output.\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# retrieve te credentials\n!cat ~/.config/earthengine/credentials","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\nuser_secret = \"credientials\" # Your user secret, defined in the add-on menu of the notebook editor\nrefresh_token = UserSecretsClient().get_secret(user_secret)\ncredentials = Credentials(\n       None,\n       refresh_token=refresh_token,\n       token_uri=ee.oauth.TOKEN_URI,\n       client_id=ee.oauth.CLIENT_ID,\n       client_secret=ee.oauth.CLIENT_SECRET,\n       scopes=ee.oauth.SCOPES)\n\n# Initialize the library.\n\nee.Initialize(credentials=credentials)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Analysing images **"},{"metadata":{"trusted":true},"cell_type":"code","source":"\n# Load Images ifrom earth engine\n\n\nroi = ee.Geometry.Polygon( [[[-67.32297404549217, 18.563112930177304],\n                            [-67.32297404549217,17.903121359128956],\n                            [-65.19437297127342, 17.903121359128956],\n                            [-65.19437297127342, 18.563112930177304]]],\n                         None, False);\n\npower_plants_table = ee.FeatureCollection(\"WRI/GPPD/power_plants\")\\\n.filterBounds(roi)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"EXPLORE THE PLANT TABLE"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get the plant data from earth engine into a dataframe\n\ndef plant_database(ft):\n    \n    \n    names=[]\n    fuel_types=[]\n    capa=[]\n    est_growth=[]\n    lats=[]\n    lons=[]\n    \n    \n    for feat in ft.getInfo()['features']:\n        \n        lat=feat['properties']['latitude']\n        long=feat['properties']['longitude']\n        fuel_type=feat['properties']['fuel1']\n        name=feat['properties']['name']\n        capacity=feat['properties']['capacitymw']\n        est_gr=feat['properties']['gwh_estimt']\n\n        names.append(name)\n        fuel_types.append(fuel_type)\n        lats.append(lat)\n        lons.append(long)\n        capa.append(capacity)\n        est_growth.append(est_gr)\n    \n    data=pd.DataFrame({'latittude':lats,'longitude':lons,'name':names,'fuel_type':fuel_types,'capacity':capa,'estimated_growth':est_growth})\n    \n    return data\n\n\n\n\nplant_datab=plant_database(power_plants_table)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plant_datab.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot the fuel type by the sum of capacity\nplant_datab.groupby(['fuel_type']).sum()['capacity'].plot(kind='bar')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Oil has the highest electricity generation , followed by Gas\nThis would mean that Oil plant factories would produced more emissions than the others"},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot the distribution of electricity generation with the plant unit location\nplant_datab.plot(kind=\"scatter\", x=\"longitude\", y=\"latittude\",\n    s=plant_datab['capacity'], label=\"electricity generated\",\n    c=\"estimated_growth\", cmap=plt.get_cmap(\"jet\"),\n    colorbar=True, alpha=0.4, figsize=(10,7),\n)\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"the Two plant units with fairely red painting has the highest electricity estimates per hour but with medium sized capacity, possibly this is a plant unit with Oil as fuel type\n\nMore than half of the plant units have about the same capacity and electricity estimates\n\nThe largest points which are in blue circles, indicate that plant units with large capacity has less electricity estimate per hour\n\nThere are points on the chart that has very small capacity and also very small electricity estimates, this may be due to the fuel type used by the plant unit"},{"metadata":{"trusted":true},"cell_type":"code","source":"plant_datab.groupby(['capacity','fuel_type'])['estimated_growth'].sum().plot(kind='barh',figsize=(10,10))\nplt.xlabel('estimated growth');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"looks like plant units with the same fuel type produces about the same Electricity\n\n**FUEL TYPE CAN BE A GOOD CHOICE TO USE WHEN ASSIGNING CLASSES TO THE FEATURE TABLE FOR CLASSIFICATION**"},{"metadata":{},"cell_type":"markdown","source":"Exploring the satelite images"},{"metadata":{"trusted":true},"cell_type":"code","source":"# code for adding layer to Map\n\n\n# Define a method for displaying Earth Engine image tiles on a folium map.\ndef add_ee_layer(self, ee_object, vis_params, name):\n    \n    try:    \n        # display ee.Image()\n        if isinstance(ee_object, ee.image.Image):    \n            map_id_dict = ee.Image(ee_object).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n            ).add_to(self)\n        # display ee.ImageCollection()\n        elif isinstance(ee_object, ee.imagecollection.ImageCollection):    \n            ee_object_new = ee_object.mosaic()\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n            ).add_to(self)\n        # display ee.Geometry()\n        elif isinstance(ee_object, ee.geometry.Geometry):    \n            folium.GeoJson(\n            data = ee_object.getInfo(),\n            name = name,\n            overlay = True,\n            control = True\n        ).add_to(self)\n        # display ee.FeatureCollection()\n        elif isinstance(ee_object, ee.featurecollection.FeatureCollection):  \n            ee_object_new = ee.Image().paint(ee_object, 0, 2)\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n        )\n        # display ee.FeatureCollection()\n        elif isinstance(ee_object, ee.feature.Feature):  \n            ee_object_new = ee.Image().paint(ee_object, 0, 2)\n            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)\n            folium.raster_layers.TileLayer(\n            tiles = map_id_dict['tile_fetcher'].url_format,\n            attr = 'Google Earth Engine',\n            name = name,\n            overlay = True,\n            control = True\n        ).add_to(self)\n    \n    except:\n        print(\"Could not display {}\".format(name))\n    \n\n\n# Add EE drawing method to folium.\nfolium.Map.add_ee_layer = add_ee_layer\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Satellite images Analysis and visualization**\n\n Map a function to generate monthly time slices per month for the 12-month period\n \n Data for Puerto Rico from July 2018 to July 2019. Exported from Earth Engine"},{"metadata":{},"cell_type":"markdown","source":"**DERVING THE NO2 EMISSIONS FROM THE SENINEL IMAGES**\n\n\nthe wind direction takes on an important role in terms of the air pollution and air quality\n\nthe tro"},{"metadata":{"trusted":true},"cell_type":"code","source":"\n# column density of NO2\ntropospheric_NO2_column_number_density=ee.ImageCollection(\"COPERNICUS/S5P/OFFL/L3_NO2\")\\\n.filterDate(\"2018-07-01\",\"2018-07-03\")\\\n.filterBounds(roi)\\\n.select('tropospheric_NO2_column_number_density')\\\n\n\n# wind data \ngfs_data=ee.ImageCollection(\"NOAA/GFS0P25\")\\\n.filterDate(\"2018-07-01\",\"2018-07-03\")\\\n.filterBounds(roi)\n\n#wind direction\nwind_direction=ee.ImageCollection(\"NASA/GLDAS/V021/NOAH/G025/T3H\")\\\n.filterDate(\"2018-07-01\",\"2018-07-03\")\\\n.filterBounds(roi)\\\n.select('Wind_f_inst')\\\n\n# temperature\ntemperature_data=ee.ImageCollection('NOAA/GFS0P25')\\\n.filterDate(\"2018-07-01\",\"2018-07-03\")\\\n.filterBounds(roi)\\\n.select('temperature_2m_above_ground')\\\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# display the two days of NO2 emission data\nmax = ee.ImageCollection(tropospheric_NO2_column_number_density).mean().reduceRegion(ee.Reducer.max(), roi, 3000).getInfo()['tropospheric_NO2_column_number_density']\n\n\nmin = ee.ImageCollection(tropospheric_NO2_column_number_density).mean().reduceRegion(ee.Reducer.min(), roi, 3000).getInfo()['tropospheric_NO2_column_number_density']\n\n\n\nvizAnomaly = {\n'min':min, 'max':max, \n'palette': ','.join(['#ffffcc','#ffeda0','#fed976','#feb24c','#fd8d3c','#fc4e2a','#e31a1c','#bd0026','#800026'])\n }\n\n\n\n# Create a folium map object.\nlat=18.200178; lon=-66.664513\nmy_map = folium.Map(location=[lat, lon],tiles=\"OpenStreetMap\"  ,  zoom_start=9.4, height=700)\nmy_map.add_ee_layer(tropospheric_NO2_column_number_density.mean().clip(roi),vizAnomaly,'NO2 emission')\nfolium.TileLayer(opacity=0.42).add_to(my_map)\n# Display the map.\n\nfor g in plant_datab.itertuples():\n\n    folium.Marker(\n    location=[g[1],g[2]],\n    popup=str(g[3]),\n    icon=folium.Icon(color='red', icon='bolt')).add_to(my_map)\n\ndisplay(my_map)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# export the latitude, longitude and array\ndef NO2_desnity_to_array(img):\n    img = img.addBands(ee.Image.pixelLonLat())\n \n    img = img.reduceRegion(reducer=ee.Reducer.toList(),\\\n                                        geometry=roi,\\\n                                        maxPixels=1e13,\\\n                                        scale=1000);\n     \n    \n    no2 = np.array((ee.Array(img.get(\"tropospheric_NO2_column_number_density\")).getInfo())) \n       \n    lats = np.array((ee.Array(img.get(\"latitude\")).getInfo()))\n    lons = np.array((ee.Array(img.get(\"longitude\")).getInfo()))\n    \n    \n    df=pd.DataFrame({'lats':lats,'lons':lons,'no2_density':no2})\n    \n    return df\n\n\n\nNO2_density=NO2_desnity_to_array(tropospheric_NO2_column_number_density.mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n#convert data into imagge\n\n# covert the lat, lon and array into an image\ndef toImage(data,image_name):\n    img_data=data[image_name]\n    # get the unique coordinates\n    \n    lats=data['lats']\n    lons=data['lons']\n    \n    uniqueLats = np.unique(data['lats'])\n    uniqueLons = np.unique(data['lons'])\n \n    # get number of columns and rows from coordinates\n    ncols = len(uniqueLons)\n    nrows = len(uniqueLats)\n \n    # determine pixelsizes\n    ys = uniqueLats[1] - uniqueLats[0]\n    xs = uniqueLons[1] - uniqueLons[0]\n \n    # create an array with dimensions of image\n    arr = np.zeros([nrows, ncols], np.float32) #-9999\n \n    # fill the array with values\n    counter =0\n    for y in range(0,len(arr),1):\n        for x in range(0,len(arr[0]),1):\n            if lats[counter] == uniqueLats[y] and lons[counter] == uniqueLons[x] and counter < len(lats)-1:\n                counter+=1\n                arr[len(uniqueLats)-1-y,x] = img_data[counter] # we start from lower left corner\n    return arr\n ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate wind speed and wind direction form \n\n# U = u_component_of_wind  , V = v_component_of_wind \n\n#wind speed = np.sqrt(U**2+ V**2)\n\n# Calculate wind direction in radians:\n\n#wind direction = np.arctan2(V,U)\n\n\ndef wind_speed_and_direction_mapper(imageCollection):\n    # u compoment of wind\n    U = imageCollection.select('u_component_of_wind_10m_above_ground')\n    # v component of wind\n    V = imageCollection.select('v_component_of_wind_10m_above_ground')\n    \n    U_component=ee.Image(U)\n    \n    V_component=ee.Image(V)\n    \n    wind_sped= U_component.pow(2).add( V_component.pow(2) ).sqrt()\\\n                                .set('system:time_start',imageCollection.get('system:time_start') )\\\n                                .rename(\"wind_speed\")\n    \n    wind_direction= V_component.atan2(U_component)\\\n                                .set('system:time_start',imageCollection.get('system:time_start') ) \\\n                                .rename(\"wind_direction\")\n    U=U.set('system:time_start',imageCollection.get('system:time_start') )\n    V=V.set('system:time_start',imageCollection.get('system:time_start') )\n    \n    return ee.Image(1).addBands([wind_sped,wind_direction,U.rename('U') ,V.rename('V')])\n\n\n\n\nwind_data=ee.ImageCollection(gfs_data).map(wind_speed_and_direction_mapper)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get the data fro wind speed and wind direction into a dataframe\n\n\n# export the latitude, longitude and array\ndef LatLonResults(img):\n    img = img.addBands(ee.Image.pixelLonLat())\n \n    img = img.reduceRegion(reducer=ee.Reducer.toList(),\\\n                                        geometry=roi,\\\n                                        maxPixels=1e13,\\\n                                        scale=1000);\n \n    wind_speed = np.array((ee.Array(img.get(\"wind_speed\")).getInfo()))\n    wind_direction = np.array((ee.Array(img.get(\"wind_direction\")).getInfo()))\n    \n    U=np.array((ee.Array(img.get(\"U\")).getInfo()))\n    V=np.array((ee.Array(img.get(\"V\")).getInfo()))\n    \n    lats = np.array((ee.Array(img.get(\"latitude\")).getInfo()))\n    lons = np.array((ee.Array(img.get(\"longitude\")).getInfo()))\n    \n    df=pd.DataFrame({'lats':lats,'lons':lons,'wind_speed':wind_speed,'wind_direction':wind_direction,'U':U,'V':V })\n    \n    return df\n\n\n\ndata=LatLonResults(wind_data.mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plotting wind speed against no2\n\n# plot the mean wind_speeds angainst the mean no2 values\nsns.regplot(x=data['wind_speed'], y=NO2_density['no2_density'],x_jitter=.05,x_estimator=np.mean);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"the regression line shows a negative slop\n\nthis means an increase in wind speed decreases the no2 tropospher column density\n\nand this is very true as pollution levels due to gas emissions in an environment are likely to be higher on low wind speed days"},{"metadata":{"trusted":true},"cell_type":"code","source":"# plotting wind direction against no2\n\n# plot the mean wind_speeds angainst the mean no2 values\nsns.regplot(x=data['wind_direction'], y=NO2_density['no2_density'],x_jitter=.05,x_estimator=np.mean);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"the regression line is nearly level but slightly tilted negatively\n\nwhich infers that high increase in wind direction contribute to a small decrease in no2 troposphere column density\n\nwind direction has a strong correlation with no2 and hence will be used in the classification process to isolate polution due to the plant units"},{"metadata":{"trusted":true},"cell_type":"code","source":"# temperature data\n\n\n# export the latitude, longitude and array\ndef temperature_to_array(img):\n    img = img.addBands(ee.Image.pixelLonLat())\n \n    img = img.reduceRegion(reducer=ee.Reducer.toList(),\\\n                                        geometry=roi,\\\n                                        maxPixels=1e13,\\\n                                        scale=1000);\n     \n    \n    no2 = np.array((ee.Array(img.get(\"temperature_2m_above_ground\")).getInfo())) \n       \n    lats = np.array((ee.Array(img.get(\"latitude\")).getInfo()))\n    lons = np.array((ee.Array(img.get(\"longitude\")).getInfo()))\n    \n    \n    df=pd.DataFrame({'lats':lats,'lons':lons,'temerature':no2})\n    \n    return df\n\n\n\ntemp=temperature_to_array(temperature_data.mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp_img=toImage(temp,'temerature')\n\n\nfig=plt.figure(figsize=(8,6))\nplt.imshow(temp_img)\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot temperature 2 above ground against no2 concentration\n\nsns.regplot(x=temp['temerature'], y=NO2_density['no2_density'],x_jitter=.05,x_estimator=np.mean);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"the regression line shows a slightly negative slop\n\nthis means an increase in temperature decreases the no2 concentration 2m above ground\n\nand this is very true as pollution levels due to gas emissions in an environment are likely to be higher on low temperature days\n\nthe scatter of the mean no2 concentration is more towards temperatures from 26 degrees onward\n\nthis shows that no2 is produced more when temperature is higher\n\nso temperature is also a good indicator for no2 pollution \n"},{"metadata":{},"cell_type":"markdown","source":"For Temperature with Plant units"},{"metadata":{"trusted":true},"cell_type":"code","source":"# get temperature of each plant unit and plot with electricity generation \n\ntemp_data=[]\nlats=[]\nlons=[]\ncapa=[]\nfuel_type=[]\nest_growth=[]\nfor g in plant_datab.itertuples():\n\n    lats.append(g[1])\n    lons.append(g[2])\n    fuel_type.append(g[4])\n    capa.append(g[5])\n    est_growth.append(g[6])\n    location=[g[2],g[1]]\n    \n    area=ee.Geometry.Point(location)\n    \n    info=ee.ImageCollection(temperature_data).mean().reduceRegion(ee.Reducer.sum(), area, 1000).getInfo()['temperature_2m_above_ground']\n    \n    \n    temp_data.append(info)\n    \n\n    \n    \ntemperature_plant_dt=pd.DataFrame({'lats':lats,'lons':lons,'fuel_type':fuel_type,'capacity':capa,'est_growth':est_growth,'temp':temp_data})    \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temperature_plant_dt.plot(kind=\"scatter\", x=\"lons\", y=\"lats\",\n    s=temperature_plant_dt['capacity'], label=\"temperature\",\n    c=\"temp\", cmap=plt.get_cmap(\"jet\"),\n    colorbar=True, alpha=0.4, figsize=(10,7),\n)\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most of te plant units here produced very high temperature when their electricity generating capacity is also high\n\nJust a few plant units with low capacity gave a low temperature this might be due to the fact that plant units use different fuel type"},{"metadata":{"trusted":true},"cell_type":"code","source":"temperature_plant_dt.groupby(['fuel_type']).sum()['temp'].plot(kind='bar')\nplt.ylabel('temperature');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"On this selected date range, Hydro and Gas gave the cumulative highest temperatures"},{"metadata":{"trusted":true},"cell_type":"code","source":"temperature_plant_dt.groupby(['capacity','fuel_type'])['temp'].sum().plot(kind='barh',figsize=(10,10))\nplt.xlabel('Temperature');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"All plant units gave a fiarly equal temperature regardless of their capacity with one plant unit of fuel type Has Given the highest temperature "},{"metadata":{},"cell_type":"markdown","source":"*For NO2 with plant units*"},{"metadata":{"trusted":true},"cell_type":"code","source":"# get no2 concentration of each plant unit and plot with electricity generation \n\nno2_data=[]\nlats=[]\nlons=[]\ncapa=[]\nfuel_type=[]\nest_growth=[]\nfor g in plant_datab.itertuples():\n\n    lats.append(g[1])\n    lons.append(g[2])\n    fuel_type.append(g[4])\n    capa.append(g[5])\n    est_growth.append(g[6])\n    location=[g[2],g[1]]\n    \n    area=ee.Geometry.Point(location)\n    \n    info=ee.ImageCollection(tropospheric_NO2_column_number_density)\\\n                                        .mean().reduceRegion(ee.Reducer.sum(), area, 1000)\\\n                                        .getInfo()['tropospheric_NO2_column_number_density']\n    \n    \n    no2_data.append(info)\n    \n\n    \n    \nno2_data_plant_dt=pd.DataFrame({'lats':lats,'lons':lons,'fuel_type':fuel_type,'capacity':capa,'est_growth':est_growth,'no2':no2_data})    \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_data_plant_dt.plot(kind=\"scatter\", x=\"lons\", y=\"lats\",\n    s=no2_data_plant_dt['capacity'], label=\"NO2 concentration\",\n    c=\"no2\", cmap=plt.get_cmap(\"jet\"),\n    colorbar=True, alpha=0.4, figsize=(10,7),\n)\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"fairly, the plant units for the selected date of NO2 emission was low, even though their capacity for high\n\njust two of those plant units with low capacity gave very high NO2 emission"},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_data_plant_dt.groupby(['fuel_type']).sum()['no2'].plot(kind='bar')\nplt.ylabel('NO2 concentration');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Again, Hydrogen and Gas Gave the highest cumulative NO2 concentration on this day, THis may be due to the number of Hydrogen and Gas plants compared to the other plants"},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_data_plant_dt.groupby(['capacity','fuel_type'])['no2'].sum().plot(kind='barh',figsize=(10,10))\nplt.xlabel('NO2 concentration');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**TEMPERATURE , WIND SPEED, WIND DIRECTION HAS AN INFLUENCE ON NO2 CONCENTRATION**"},{"metadata":{},"cell_type":"markdown","source":"**USE A MONTHS DATA TO EXPERIMENT WITH THE METHODOLOGY**"},{"metadata":{},"cell_type":"markdown","source":" The project is to calculate the emission factor with respect to only electricity generation\n \n It is True that Puerto rico has a well distinct energy source which will be a high contribution factor to NO2 concentration\n \n But other factors are also there in Puerto Rico which contribute to NO2 Concentration such as traffic, e.t.c.\n \n Therefore we have to use a statistical model to distinguish NO2 concentration sources and isolate concentrations of NO2 attributable to electricity\n \n A statistical  method  in  analyzing  the  air pollution  sources  has  been  widely  used  using  the  Principal  Component  Analysis  (PCA) \n \n**PCA is a one of the multivariate statistical analysis which categorizing the data into similar group, provided it has relationship of each other and later known as a principal component.  **"},{"metadata":{"trusted":true},"cell_type":"code","source":"# wind data \ngfs_data=ee.ImageCollection(\"NOAA/GFS0P25\")\\\n.filterDate(\"2018-07-25\",\"2018-07-30\")\\\n.filterBounds(roi)\n\n\nwind_data=ee.ImageCollection(gfs_data).map(wind_speed_and_direction_mapper)\n\n# temperature\ntemp_d=ee.ImageCollection('NOAA/GFS0P25')\\\n.filterDate(\"2018-07-25\",\"2018-07-30\")\\\n.filterBounds(roi)\\\n.select('temperature_2m_above_ground')\\\n\n\n# NO2 imagecollection\nNO2=ee.ImageCollection(\"COPERNICUS/S5P/OFFL/L3_NO2\")\\\n.filterDate(\"2018-07-01\",\"2018-07-30\")\\\n.filterBounds(roi)\\\n.select('tropospheric_NO2_column_number_density')\\\n\nselector=['tropospheric_NO2_column_number_density','temperature_2m_above_ground','wind_speed','wind_direction']\n\n\n\nninput=NO2.mean().addBands([temp_d.mean(),wind_data.mean()]).select(selector)\n\n\n#************************************************* PCA************************************************\n\n# Get some information about the input to be used later.\nscale = ninput.projection().nominalScale().getInfo()\n\nbandNames = ninput.bandNames()\n\n#Mean center the data to enable a faster covariance reducer and an SD stretch of the principal components.\n#reducer, geometry, scale, crs, crsTransform, bestEffort, maxPixels, tileScale\nmeanDict = ninput.reduceRegion(ee.Reducer.mean(),roi,scale,maxPixels=1e9)\n\nmeans = ee.Image.constant(meanDict.values(bandNames))\ncentered = ninput.subtract(means)\n\n\n\n\n\n#This helper function returns a list of new band names.\ndef getNewBandNames(prefix):\n    seq = ee.List.sequence(1, bandNames.length())    \n    mlist=[]\n    for x in seq.getInfo():\n        new=f\"{prefix}{x}\"\n        mlist.append(new)\n    return ee.List(mlist)\n                 \n# This function accepts mean centered imagery, a scale and a region in which to perform the analysis.  \n# It returns the Principal Components (PC) in the region as a new image.\n                  \ndef getPrincipalComponents(centered, scale, region):\n    #Collapse the bands of the image into a 1D array per pixel.\n    arrays = centered.toArray()\n  \n    #Compute the covariance of the bands within the region.\n    covar = arrays.reduceRegion(ee.Reducer.centeredCovariance(),region, scale,maxPixels=1e9)\n    \n    #Get the 'array' covariance result and cast to an array. This represents the band-to-band covariance within the region.\n    covarArray = ee.Array(covar.get('array'))\n  \n    #Perform an eigen analysis and slice apart the values and vectors.\n    eigens = covarArray.eigen()\n    \n    #This is a P-length vector of Eigenvalues.\n    eigenValues = eigens.slice(1, 0, 1)\n    \n    #This is a PxP matrix with eigenvectors in rows.\n    eigenVectors = eigens.slice(1, 1)\n    \n    #Convert the array image to 2D arrays for matrix computations.\n    arrayImage = arrays.toArray(1)\n    \n    #Left multiply the image array by the matrix of eigenvectors.\n    principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage)\n    \n    \n    #Turn the square roots of the Eigenvalues into a P-band image.\n    sdImage = ee.Image(eigenValues.sqrt()).arrayProject([0]).arrayFlatten([getNewBandNames('sd')])\n    #Turn the PCs into a P-band image, normalized by SD.\n    \n    return principalComponents.arrayProject([0]).arrayFlatten([getNewBandNames('pc')]).divide(sdImage)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Get the PCs at the specified scale and in the specified region\npcImage = getPrincipalComponents(centered, scale, roi)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"band = pcImage.bandNames().get(0).getInfo()\nprint(band)\nprint(pcImage.select([band]).reduceRegion(ee.Reducer.minMax(),roi,1000).getInfo())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"the PCA images for NO2 concentration, temperature, wind speed and wind direction gives positive and neative values of PCA\n\nWhen different variables are subjected to PCA, not all are positively correlated with each other. There are some variables which are negatively correlated. That's why, theoretically variance as represented by e-values in PCA may b negative or positive\n\nPCA 1 = variance of other variables with respect to NO2 concentration\n\nPCA 2 = variance with of other variables with respect to temperature\n\nPCA 3 = variance with of other variables with respect to wind speed\n\nPCA 4 = variance with of other variables with respect to wind direction\n\n\nTherefore for PCA 1,2,3,4 , positive values represent areas that are positively correlated in terms of NO2 concentration, temperature, wind speed, and wind direction.\n\nThis gives us a fair idea of pollution sources\n\nhence  highly positive correlated areas can represent pollution due to the power plants\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"#*****************   PCA1  tropospheric_NO2_column_number_density\n \n# selector=['tropospheric_NO2_column_number_density','temperature_2m_above_ground','wind_speed','wind_direction']\n\n# band = pcImage.bandNames().get(0).getInfo()\n# print(band)\n# print(pcImage.select([band]).reduceRegion(ee.Reducer.minMax(),roi,1000))\n \n#Map.addLayer(pcImage.select([band]), {min:-3.2, max:2 ,palette:p1}, band);\n\n\nvizAnomaly = {\n'min':-1.3491126012357098, 'max': 6.6613191472965845, \n'palette': ','.join(['#ffffcc','#ffeda0','#fed976','#feb24c','#fd8d3c','#fc4e2a','#e31a1c','#bd0026','#800026'])\n }\n\n\n\n# Create a folium map object.\nlat=18.200178; lon=-66.664513\nmy_map = folium.Map(location=[lat, lon],tiles=\"OpenStreetMap\"  ,  zoom_start=9.4, height=700)\nmy_map.add_ee_layer(pcImage.select(['pc1']).clip(roi),vizAnomaly,'pc1')\nfolium.TileLayer(opacity=0.1).add_to(my_map)\n\n#add the power plants location\nfor g in plant_datab.itertuples():\n\n    folium.Marker(\n    location=[g[1],g[2]],\n    popup=str(g[3]),\n    icon=folium.Icon(color='blue', icon='info-sign')).add_to(my_map)\n\n\n\n# Display the map.\ndisplay(my_map)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From PCA 1 image, we see that places with not power plant have very high PCA than places wit power plants, \n\nthis might be due to pollution due to population and traffic or other pollution sources other than the power plants electricity generation"},{"metadata":{"trusted":true},"cell_type":"code","source":"band = pcImage.bandNames().get(1).getInfo()\nprint(band)\nprint(pcImage.select([band]).reduceRegion(ee.Reducer.minMax(),roi,1000).getInfo())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#*****************PCA2 temperature\n\n# band = pcImage.bandNames().get(1).getInfo()\n \n# print(pcImage.select([band]).reduceRegion(ee.Reducer.minMax(),roi,1000))\n\n\nvizAnomaly = {\n'min':-2.8838120828363216, 'max': 24.497258166422363, \n'palette': ','.join(['#ffffcc','#ffeda0','#fed976','#feb24c','#fd8d3c','#fc4e2a','#e31a1c','#bd0026','#800026'])\n }\n\n\n\n# Create a folium map object.\nlat=18.200178; lon=-66.664513\nmy_map = folium.Map(location=[lat, lon],tiles=\"OpenStreetMap\"  ,  zoom_start=9.4, height=700)\nmy_map.add_ee_layer(pcImage.select(['pc2']).clip(roi),vizAnomaly,'pc2')\nfolium.TileLayer(opacity=0.1).add_to(my_map)\n# Display the map.\n\n#add the power plants location\nfor g in plant_datab.itertuples():\n\n    folium.Marker(\n    location=[g[1],g[2]],\n    popup=str(g[3]),\n    icon=folium.Icon(color='red', icon='info-sign')).add_to(my_map)\n\n\n\ndisplay(my_map)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From PCA 2 image, we see that places with no or less number of power plant have low  PCA than places with power plants, \n\nThe temperatures from power plants are higher than that of other pollution sources like traffic e.t.c\n\nThe temperatures might also differ from power plants due to the different fuel types\n\n**hence PCA 2 is a good predictor of pollution source**\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"band = pcImage.bandNames().get(2).getInfo()\nprint(band)\nprint(pcImage.select([band]).reduceRegion(ee.Reducer.minMax(),roi,1000).getInfo())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#*****************PCA3\n\nvizAnomaly = {\n'min': -215968294.192117, 'max': 32334311.396379523, \n'palette': ','.join(['#ffffcc','#ffeda0','#fed976','#feb24c','#fd8d3c','#fc4e2a','#e31a1c','#bd0026','#800026'])\n }\n\n\n \n# Create a folium map object.\nlat=18.200178; lon=-66.664513\nmy_map = folium.Map(location=[lat, lon],tiles=\"OpenStreetMap\"  ,  zoom_start=9.4, height=700)\nmy_map.add_ee_layer(pcImage.select(['pc3']).clip(roi),vizAnomaly,'pc3')\nfolium.TileLayer(opacity=0.1).add_to(my_map)\n\n#add the power plants location\nfor g in plant_datab.itertuples():\n\n    folium.Marker(\n    location=[g[1],g[2]],\n    popup=str(g[3]),\n    icon=folium.Icon(color='green', icon='info-sign')).add_to(my_map)\n\n\n# Display the map.\ndisplay(my_map)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From PCA 3 image, we see that the distribution of PCA scores is not really given us the information we need to know pollution sources\n\nBecause even areas which no power plant can either have high or low PCA score, **hence PCA 3 is not a good predictor of pollution source**"},{"metadata":{},"cell_type":"markdown","source":"USING PCA 2 IMAGE TO MUSK OUT PIXELX IN NO2 TROPOSPHERE VERTICAL COLUMN TO GET NO2 CONCENTRATIONS ATTRIBUTABLE TO POWER PLANTS"},{"metadata":{"trusted":true},"cell_type":"code","source":"# pca 2 image\npca_2=pcImage.select(['pc2'])\n\n#no2 concentration image\nno2_img=NO2.mean()\n\n# since pca2 gave both positive and negative pca scores, we want to musk out pixels in no2 image that had negative variation in pca2 image\n\n#Create a binary mask.\nmask = pca_2.gt(0) # pixel greater than zero\n\n# Update the composite mask with the pca 2 mask.\nmaskedComposite = no2_img.updateMask(mask)\n\nprint(maskedComposite.getInfo())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"max = maskedComposite.reduceRegion(ee.Reducer.max(), roi, 3000).getInfo()['tropospheric_NO2_column_number_density']\n\n\nmin = maskedComposite.reduceRegion(ee.Reducer.min(), roi, 3000).getInfo()['tropospheric_NO2_column_number_density']\n\n\nprint(min,max)\n\nvizAnomaly = {\n'min':min, 'max':max, \n'palette': ','.join(['#e7e1ef','#c994c7','#dd1c77'])\n }\n\n\n\n# Create a folium map object.\nlat=18.200178; lon=-66.664513\nmy_map = folium.Map(location=[lat, lon],tiles=\"OpenStreetMap\"  ,  zoom_start=9.4, height=700)\nmy_map.add_ee_layer(maskedComposite.clip(roi),vizAnomaly,'NO2 concentration')\nfolium.TileLayer(opacity=0.1).add_to(my_map)\n# Display the map.\n\nfor g in plant_datab.itertuples():\n\n    folium.Marker(\n    location=[g[1],g[2]],\n    popup=str(g[3]),\n    icon=folium.Icon(color='red', icon='info-sign')).add_to(my_map)\n\ndisplay(my_map)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The image displayed gives a fiar image of no2 concentrations due to the power plants"},{"metadata":{},"cell_type":"markdown","source":"NO2 emission = NO2 Concentration * Area occupied by the NO2 concentrations"},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_concentration = maskedComposite.clip(roi)\n\n# tropospheric no2 vertical column in mol/m^2 ->\n# calculate the pixel area in square meters\narea_no2 = no2_concentration.multiply(ee.Image.pixelArea()).divide(1000*1000).multiply(1000000)\n\n\n\n\n# reducing the statistics for the area\nstat = area_no2.reduceRegion(ee.Reducer.sum(),roi,1000, maxPixels=1e9).getInfo()['tropospheric_NO2_column_number_density']\n\n#get the sq meter area for no2 concentration\nprint ('NO2 concentration Area (in sq.meter) is', stat)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_emission= no2_concentration.multiply(stat)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"max = no2_emission.reduceRegion(ee.Reducer.max(), roi, 3000).getInfo()['tropospheric_NO2_column_number_density']\n\n\nmin = no2_emission.reduceRegion(ee.Reducer.min(), roi, 3000).getInfo()['tropospheric_NO2_column_number_density']\n\n\nprint(min,max)\n\nvizAnomaly = {\n'min':min, 'max':max, \n'palette': ','.join(['#ffeda0','#feb24c','#f03b20'])\n }\n\n\n\n# Create a folium map object.\nlat=18.200178; lon=-66.664513\nmy_map = folium.Map(location=[lat, lon],tiles=\"OpenStreetMap\"  ,  zoom_start=9.4, height=700)\nmy_map.add_ee_layer(no2_emission.clip(roi),vizAnomaly,'NO2 concentration')\nfolium.TileLayer(opacity=0.1).add_to(my_map)\n# Display the map.\n\nfor g in plant_datab.itertuples():\n\n    folium.Marker(\n    location=[g[1],g[2]],\n    popup=str(g[3]),\n    icon=folium.Icon(color='red', icon='info-sign')).add_to(my_map)\n\ndisplay(my_map)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**EMISSION FACOTR (EF) = EMISSION / ACTICITY **\n\nTHE ACTIVITY GENERATING THE NO2 EMISSION IS THE ELECTRICITY GENERATION FROM THE POWER PLANTS"},{"metadata":{"trusted":true},"cell_type":"code","source":"# calculate the total electricity generated\nsum_capacity=plant_datab['capacity'].sum()\n\n\nelectricity_gen= sum_capacity","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Emission factor\n\nemission_factor= no2_emission.divide(electricity_gen)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# reducing the statistics for the area\nstat = emission_factor.reduceRegion (ee.Reducer.sum(),roi,1000, maxPixels=1e9).getInfo()['tropospheric_NO2_column_number_density']\n\n\nprint ('Emission factor for sub_urban area of peurto rico is', stat)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**MARGINAL EMISSION FACTOR**\n\nThe correct way to measure the impact of environmental decisions is to use marginal emissions factors.\n\nMarginal emissions factors measure the actual environmental consequences of taking different potential actions on the power grid.\n\nIf the example city is evaluating an energy efficiency measure to conserve one megawatt-hour ofelectricity consumption, this program will reduce local emissions by reducing output at one or more power plants. But which power plants? Many sources of power, for example most solar\n\n**Conserving energy only affects some power plants: those which can scale up or down in response, known as the “marginal” power plants.**\n\n**Marginal emissions measure the emissions per kilowatt-hour only from these power plants, thus accurately measuring real-world results.**\n\n\n\n**CALCULATING MARGINAL EMISSION FACTOR**\n\n   * FIRST CALCULATE THE EFFICIENCY OF ELECTRICITY GENERATION POR EACH POWER UNIT*\n    \n   * SUM ALL THE EFFICIENCY OF THE POWER UNITS TO GET THE NET EFFICIENCY*\n    \n   * DIVIDE THE MEAN NO2 EMISION BY THE EFFIEINCY*\n    \n**MEF = MEAN NO2 EMISSION / NET EFFICIENCY  **  "},{"metadata":{"trusted":true},"cell_type":"code","source":"#CALCULATE THE EFFICIENCY OF ELECTRICITY GENERATION POR EACH POWER UNIT\n# EFFICIENCY = ELECTRICITY CAPACITY - ESTIMATED ELECTRICITY\n\n# convert megawatt hour to gigawatt hours \ncapacity_elect=(plant_datab['capacity'] /1000) * 24 \n\nestimated_elect= plant_datab['estimated_growth']\n\nplant_datab['efficiency']= estimated_elect - capacity_elect \n\nnet_efficiency = plant_datab.efficiency.sum()\n\nprint ('Net Efficiency for power plants of peurto rico is', f'{net_efficiency} gigawatt hours')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# CALCULATING MARGINAL EMISSION FACTOR\n\nemf= no2_emission.divide(net_efficiency)\n\n# reducing the statistics for the area\nstat = emf.reduceRegion (ee.Reducer.sum(),roi,1000, maxPixels=1e9).getInfo()['tropospheric_NO2_column_number_density']\n\n\nprint ('Marginal Emission factor for sub_urban area of peurto rico is', stat)\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**NOW CALCULATE THE AVERAGE HISTORICAL EMISSION FACTOR AND MARGINAL EMISSION FACTOR USING THE ABOVE PROCEDURES FOR PUERTO RICO**\n\nthis time get the monthly 3 days data for temperature and wind data and average them and recompute the emision factor, this is because the gfs data is huge and always gives memory issuse"},{"metadata":{"trusted":true},"cell_type":"code","source":"\n#calculate monthly data for temperature\n# because of memory issues from google earth engine api, get the last 10 days data for each month and average it\ndef calcMonthlyMean_temp(years,months):\n    mylist = ee.List([])\n    for y in years:\n        for m in months:\n            if(y==2018 and m>6 ) or (y==2019 and m<8) :\n                \n                \n                gfs_imagecollection=ee.ImageCollection(\"NOAA/GFS0P25\")\\\n                                .filter(ee.Filter.date(ee.Date.fromYMD(y,m,1),ee.Date.fromYMD(y,m,3)))\\\n                                .filterBounds(roi)\\\n                                .select('temperature_2m_above_ground')\n              \n                \n                w = gfs_imagecollection.filter(ee.Filter.calendarRange(y, y, 'year'))\\\n                .filter(ee.Filter.calendarRange(m, m, 'month')).mean()\n            \n                mylist = mylist.add(w.set('year', y)\\\n                                    .set('month', m)\\\n                                    .set('date', ee.Date.fromYMD(y,m,1))\\\n                                    .set('system:time_start',\\\n                                         ee.Date.fromYMD(y,m,1)))\n            \n    return ee.ImageCollection.fromImages(mylist)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"years = [2018, 2019]\n\nmonths = range(1,13)\n\n# get monthly temperature data from gfs data\nmonthly_temperature= calcMonthlyMean_temp(years,months)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#calculate monthly data for wind speed and wind direction \n# because of memory issues from google earth engine api, get the last 10 days data for each month and average it\ndef calcMonthlyMean_wind(years,months):\n    mylist = ee.List([])\n    for y in years:\n        for m in months:\n            if(y==2018 and m>6 ) or (y==2019 and m<8) :\n                \n                \n                gfs_imagecollection=ee.ImageCollection(\"NOAA/GFS0P25\")\\\n                                .filter(ee.Filter.date(ee.Date.fromYMD(y,m,1),ee.Date.fromYMD(y,m,3)))\\\n                                .filterBounds(roi)\\\n                                .select(['u_component_of_wind_10m_above_ground','v_component_of_wind_10m_above_ground'])\n\n              \n                \n                w = gfs_imagecollection.filter(ee.Filter.calendarRange(y, y, 'year'))\\\n                .filter(ee.Filter.calendarRange(m, m, 'month')).mean()\n            \n                mylist = mylist.add(w.set('year', y)\\\n                                    .set('month', m)\\\n                                    .set('date', ee.Date.fromYMD(y,m,1))\\\n                                    .set('system:time_start',\\\n                                         ee.Date.fromYMD(y,m,1)))\n            \n    return ee.ImageCollection.fromImages(mylist)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# wind data \n\nyears = [2018, 2019]\n\nmonths = range(1,13)\n\n# get monthly wind speed and wind direction data from gfs data\nmonthly_wind_s_d= calcMonthlyMean_wind(years,months)\n                ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#no2 concentration\nno2=ee.ImageCollection(\"COPERNICUS/S5P/OFFL/L3_NO2\")\\\n.filterDate(\"2018-07-01\",\"2019-07-01\")\\\n.filterBounds(roi)\\\n.select('tropospheric_NO2_column_number_density')\\\n\n#temperature\ntemp = monthly_temperature\n\n#wind_data\nwind = monthly_wind_s_d\nwind_data=ee.ImageCollection(wind).map(wind_speed_and_direction_mapper).select(['wind_speed','wind_direction'])\n\n\nselector=['tropospheric_NO2_column_number_density','temperature_2m_above_ground','wind_speed','wind_direction']\n\nninput=no2.mean().addBands([temp.mean(),wind_data.mean()]).select(selector)\n\n\n#************************************************* PCA ************************************************\n\n# Get some information about the input to be used later.\nscale = ninput.projection().nominalScale().getInfo()\n\nbandNames = ninput.bandNames()\n\n#Mean center the data to enable a faster covariance reducer and an SD stretch of the principal components.\n#reducer, geometry, scale, crs, crsTransform, bestEffort, maxPixels, tileScale\nmeanDict = ninput.reduceRegion(ee.Reducer.mean(),roi,scale,maxPixels=1e9)\n\nmeans = ee.Image.constant(meanDict.values(bandNames))\ncentered = ninput.subtract(means)\n\n#Get the PCs at the specified scale and in the specified region\npcImage = getPrincipalComponents(centered, scale, roi)\n\n# pca 2 image\npca_2=pcImage.select(['pc2'])\n\n#************************************************* NO2 EMISSION ************************************************                             \n#no2 concentration image\nno2_img=no2.mean()\n\n#Create a binary mask.\nmask = pca_2.gt(0) # pixel greater than zero\n\n# Update the composite mask with the pca 2 mask.\nmaskedComposite = no2_img.updateMask(mask)\n\n\nno2_concentration = maskedComposite.clip(roi)\n\n\n# tropospheric no2 vertical column in mol/m^2 ->\n# calculate the pixel area in square meters\narea_no2 = no2_concentration.multiply(ee.Image.pixelArea()).divide(1000*1000).multiply(1000000)\n\n# reducing the statistics for the area\nstat = ee.Image(area_no2).reduceRegion(ee.Reducer.sum(),roi,10000, maxPixels=1e9).getInfo()['tropospheric_NO2_column_number_density']\n\n\nno2_emission= no2.mean().multiply(stat)\n\n#************************************************* NO2 EMISSION FACTOR ************************************************                             \n#  # convert megawatt  to gigawatt hours and calculate the total electricity generated\nsum_capacity=((plant_datab['capacity'] /1000) * 24 ).sum()\n\n\nelectricity_gen= sum_capacity\n\n#Emission factor\nemission_factor= no2_emission.divide(electricity_gen).set('system:time_start',no2.get('system:time_start') ).rename(\"ef\")\n\nef_stat = emission_factor.reduceRegion (ee.Reducer.sum(),roi,10000, maxPixels=1e9).getInfo()['ef']\n\n#************************************************* NO2 MARGINAL EMISSION FACTOR ************************************************                             \n#CALCULATE THE EFFICIENCY OF ELECTRICITY GENERATION POR EACH POWER UNIT\n# EFFICIENCY = ELECTRICITY CAPACITY - ESTIMATED ELECTRICITY\n\n# convert megawatt  to gigawatt hours \ncapacity_elect=(plant_datab['capacity'] /1000) * 24 \n\nestimated_elect= plant_datab['estimated_growth']\n\nplant_datab['efficiency']= estimated_elect - capacity_elect \n\nnet_efficiency = plant_datab.efficiency.sum()\n\n# CALCULATING MARGINAL EMISSION FACTOR\n\nmef= no2_emission.divide(net_efficiency).set('system:time_start',no2.get('system:time_start') ).rename(\"mef\")\n\nmef_stat = mef.reduceRegion (ee.Reducer.sum(),roi,10000, maxPixels=1e9).getInfo()['mef']\n\n#************************************************* RESULTS ************************************************                             \n\nprint ('Average historical Emission factor for sub_urban area of peurto rico is', ef_stat)\nprint('-'*40)\nprint ('Average historical Marginal Emission factor for sub_urban area of peurto rico is', mef_stat)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**MODEL SCALABILITY AND ACCURACY**\n\nbecause all computations where done using earth engine algorithms,which enables us to process large copous of data and get instant results, this method is very scalable as it can take very large geograhpic bounding box\n\nthe prinicpal component analysis used is very accurate in identifying pollution sources in the image data and hence this method is quite accurate\n\n\n**MODEL DEPENDANCY**\n\n1.this model improves the current emission factor calculation as it is able to clearly distinguish or isolate pollution related to electricity \n\n2.No need for data collection on fuel consumption or other difficult data parameters for computing emission factor\n\n\n**RECOMMENDATION**\n\nit is recommended that to make this model more roboust, we need to use classification algorithms and more factors such as aerosel data,ozone data and other weather data to make or improve the model"}],"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}