{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Goal\n\nDevelop a methodology to calculate an **average historical emissions factor of electricity generated** for a sub-national region, using remote sensing data and techniques.  \n(from competition overview)\n\n## Emission factor\n\nAn 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. 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). 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).  \nThe general equation for emissions estimation is as follows:\n\n**E = A × EF × (1–ER)/100**  \n\nwhere E = emissions; A = activity rate; EF = emission factor, and ER = overall emission reduction efficiency, %.  \n[Source](https://www.sciencedirect.com/topics/engineering/emission-factor)"},{"metadata":{},"cell_type":"markdown","source":"**Imports**"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"!pip install rasterstats --quiet","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"## Importing necessary libraries\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\nsns.set(style=\"white\", palette=\"muted\", color_codes=True)\n\nfrom ast import literal_eval\nfrom datetime import datetime, timedelta\nimport os, random\n\n# Plotting geographical data\nimport folium\nfrom folium import plugins\nimport rasterio as rio\nimport rasterstats\nimport geopandas as gpd\nfrom shapely.geometry import Point","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Explore data\n\nTake a look at data provided."},{"metadata":{},"cell_type":"markdown","source":"## Global powerplant database\n\nStart with provided csv of power plant in Puerto Rico  \n(Note that's also available through [Google Earth Engine](https://developers.google.com/earth-engine/datasets/catalog/WRI_GPPD_power_plants). Maybe for future use?)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"global_power_plants = pd.read_csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\nglobal_power_plants.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"global_power_plants.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"global_power_plants.describe().T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some features seem to be empty (generation_gwh_year, etc.) and .geo column needs some work to be useful."},{"metadata":{},"cell_type":"markdown","source":"## Clean gpp csv\nRemove useless features"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Columns with only 0 or NaN values\nto_drop = [\"generation_gwh_2013\", \n           \"generation_gwh_2014\", \n           \"generation_gwh_2015\", \n           \"generation_gwh_2016\",\n           \"generation_gwh_2017\", \n           \"other_fuel1\",\n           \"other_fuel2\",\n           \"other_fuel3\",\n           \"year_of_capacity_data\"]\n\nglobal_power_plants = global_power_plants.drop(to_drop, axis=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As power plants are all supposed to be in Puerto Rico, country, country_long, and geolocation_source cols might be useless.  \nLet's check and remove if necessary."},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants['country'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants['country_long'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants['geolocation_source'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Columns with all same values\nto_drop = [\"country\", \n           \"country_long\", \n           \"geolocation_source\"\n          ]\n\nglobal_power_plants = global_power_plants.drop(to_drop, axis=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's get lat/lon from .geo column.  \nThere are several ways to do it. Here, let simply convert the .geo info to a dict and recover lat and lon"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def get_lon_from_geo(str_):\n    dict_ = literal_eval(str_)\n    coordinates = dict_['coordinates']\n    lon = coordinates[0]\n    return lon\n\ndef get_lat_from_geo(str_):\n    dict_ = literal_eval(str_)\n    coordinates = dict_['coordinates']\n    lat = coordinates[1]\n    return lat\n\nglobal_power_plants['lon'] = global_power_plants['.geo'].map(get_lon_from_geo)\nglobal_power_plants['lat'] = global_power_plants['.geo'].map(get_lat_from_geo)\n\nglobal_power_plants.drop(columns=['.geo'], inplace=True)\n\nglobal_power_plants.head(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's also create a GeoPandas version for future use  \n(inspired from https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson )"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"geometry_power_plants = [Point(x,y) for x,y in zip(global_power_plants['lon'], global_power_plants['lat'])]\nglobal_power_plants_gdf = gpd.GeoDataFrame(global_power_plants, crs = {'init': 'epsg: 4326'}, geometry = geometry_power_plants)\n\nglobal_power_plants_gdf.plot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Look at power plant data"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2,figsize=(12,5))\n\nsns.countplot(x=\"primary_fuel\",\n              data=global_power_plants,\n              ax=ax[0]\n             )\nsns.stripplot(x=\"primary_fuel\",\n            y=\"capacity_mw\",\n            data=global_power_plants,\n            ax=ax[1]\n           )\n\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lots of renewables, but with low capacity.  \nOil and gas power plants are of multiple capacity.\nOnly one coal power plant.\n\nLet's sum to see which technology has the most installed capacity."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sns.barplot(x='primary_fuel',\n            y=\"capacity_mw\",\n            data=global_power_plants.groupby(['primary_fuel']).sum().reset_index())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As expected, oil is clearly ahead.\n\nLet's check on generation side."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sns.barplot(x=\"primary_fuel\",\n            y=\"estimated_generation_gwh\",\n            data=global_power_plants.groupby(['primary_fuel']).sum().reset_index())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Coal is clearly out of range. Problem in dataset?  \nLet's focus on other sources."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sns.barplot(x=\"primary_fuel\",\n            y=\"estimated_generation_gwh\",\n            data=global_power_plants[global_power_plants['primary_fuel'] != 'Coal'].groupby(['primary_fuel']).sum().reset_index())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Solar and wind seem to do not contribute to energy generation.  \nLet's look in the details."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"solar_wind_power_plants = global_power_plants[(global_power_plants['primary_fuel'] == 'Solar') |\n                                              (global_power_plants['primary_fuel'] == 'Wind')]\n\nsolar_wind_power_plants","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Solar and wind seem to be neglible in term of reported generated power in the dataset.  \nWe need to see to which period correspond the value in \"estimated_generation_gwh\".\n\nAnyway, for EF, we're more interested by oil, gas and coal."},{"metadata":{},"cell_type":"markdown","source":"### Map of power plant in Puerto Rico\n\nLet's vizualize the power plants on a map."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def initMap(df, lat, lon):\n    location = [lat, lon]\n    Map = folium.Map(location=location, zoom_start=9)\n    \n    fuelColor = {'Coal': 'darkred',\n                 'Oil': 'black',\n                 'Gas': 'lightgray',\n                 'Hydro': 'lightblue',\n                 'Solar': 'orange',\n                 'Wind': 'green'\n                }\n    \n    for _, row in df.iterrows():\n        fuel = row['primary_fuel']\n        capacity = row['capacity_mw']\n        generation = row['estimated_generation_gwh']\n        locationpp = [row['lat'], row['lon']]\n        popup = \"{} {}MW capacity, {}GWh generated\".format(fuel,capacity,generation)\n        color = fuelColor[fuel]\n        folium.Marker(locationpp,\n                      popup=popup,\n                      icon=folium.Icon(color=color, icon_color='white', icon='bolt', prefix='fa')\n                     ).add_to(Map)  \n\n    return Map","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"lat=18.200178; lon=-66.664513 #puerto rico\nMap = initMap(global_power_plants, lat, lon)   \nMap","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sentinel-5P OFFL NO2: Offline Nitrogen Dioxide\n\nLet's now look at the remote sensing data, starting from Sentinel.\n\n*This dataset provides offline high-resolution imagery of NO2 concentrations.  \nNitrogen oxides (NO2 and NO) are important trace gases in the Earth’s atmosphere, present in both the troposphere and the stratosphere. They enter the atmosphere as a result of anthropogenic activities (**notably fossil fuel combustion and biomass burning**) and natural processes (w**ildfires, lightning, and microbiological processes in soils**). Here, NO2 is used to represent concentrations of collective nitrogen oxides because during daytime, i.e. in the presence of sunlight, a photochemical cycle involving ozone (O3) converts NO into NO2 and vice versa on a timescale of minutes. The TROPOMI NO2 processing system is based on the algorithm developments for the DOMINO-2 product and for the EU QA4ECV NO2 reprocessed dataset for OMI, and has been adapted for TROPOMI. This retrieval-assimilation-modelling system uses the 3-dimensional global TM5-MP chemistry transport model at a resolution of 1x1 degree as an essential element. * \n(From GEE doc)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sentinel_path = \"../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/\"\nexamples = [random.choice(os.listdir(sentinel_path)) for _ in range(4)]\nexamples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Images seem to be from time plage.\nLet's see, what they contain."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"image_name = random.choice(os.listdir(sentinel_path))\nimage_path = sentinel_path + image_name\nimage = rio.open(image_path)\n\nbands = image.count\nprint(f\"There are {bands} bands\")\n\nfor i in image.indexes:\n    desc = image.descriptions[i-1]\n    print(f'{i}: {desc}')\n\nprint('\\n')\nprint(f'Example of:{image_name}')\n\nfig, ax = plt.subplots(3,4,figsize=(20,6))\nband=1\nfor n in range(3):\n    for m in range(4):\n        image_bandi = image.read(band)\n        desc = image.descriptions[band-1]\n        ax[n,m].set_title(desc)\n        ax[n,m].imshow(image_bandi, cmap=\"Reds\")\n        ax[n,m].grid(False)\n        band += 1\n        \nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"tropospheric_NO2_column_number_density might be particulary interesting as it's related to emissions close to the ground."},{"metadata":{},"cell_type":"markdown","source":"## GFS: Global Forecast System 384-Hour Predicted Atmosphere Data\n\nLet's now look at gfs data.\n\n*The Global Forecast System (GFS) is a weather forecast model produced by the National Centers for Environmental Prediction (NCEP). The **GFS dataset consists of selected model outputs** (described below) as gridded forecast variables. The 384-hour forecasts, with 3-hour forecast interval, are made at 6-hour temporal resolution (i.e. u**pdated four times daily**). Use the 'creation_time' and 'forecast_time' properties to select data of interest.  \nThe GFS is a coupled model, composed of an atmosphere model, an ocean model, a land/soil model, and a sea ice model which work together to provide an accurate picture of weather conditions.  *  \n(From GEE)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gfs_path = \"../input/ds4g-environmental-insights-explorer/eie_data/gfs/\"\nexamples = [random.choice(os.listdir(gfs_path)) for _ in range(4)]\nexamples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Several files per day, as expected from doc.  \nLet's see how they look like."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"image_name = random.choice(os.listdir(gfs_path))\nimage_path = gfs_path + image_name\nimage = rio.open(image_path)\n\nbands = image.count\nprint(f\"There are {bands} bands\")\n\nfor i in image.indexes:\n    desc = image.descriptions[i-1]\n    print(f'{i}: {desc}')\n\nprint('\\n')\nprint(f'Example of:{image_name}')\n\nfig, ax = plt.subplots(3,2,figsize=(15,6))\nband=1\nfor n in range(3):\n    for m in range(2):\n        image_bandi = image.read(band)\n        desc = image.descriptions[band-1]\n        ax[n,m].set_title(desc)\n        ax[n,m].imshow(image_bandi, cmap=\"Reds\")\n        ax[n,m].grid(False)\n        band += 1\n        \nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## GLDAS-2.1: Global Land Data Assimilation System\n\nLet's now look at GLDAS data.\n\n*Global Land Data Assimilation System (GLDAS) ingests **satellite and ground-based observational data products**. Using advanced land surface modeling and data assimilation techniques, it generates **optimal fields of land surface states and fluxes**.*\n\n(From GEE)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gldas_path = \"../input/ds4g-environmental-insights-explorer/eie_data/gldas/\"\nexamples = [random.choice(os.listdir(gldas_path)) for _ in range(4)]\nexamples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems to have several files per day.  \nLet's vizualize some."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"image_name = random.choice(os.listdir(gldas_path))\nimage_path = gldas_path + image_name\nimage = rio.open(image_path)\n\nbands = image.count\nprint(f\"There are {bands} bands\")\n\nfor i in image.indexes:\n    desc = image.descriptions[i-1]\n    print(f'{i}: {desc}')\n\nprint('\\n')\nprint(f'Example of:{image_name}')\n\nfig, ax = plt.subplots(3,4,figsize=(20,6))\nband=1\nfor n in range(3):\n    for m in range(4):\n        image_bandi = image.read(band)\n        desc = image.descriptions[band-1]\n        ax[n,m].set_title(desc)\n        ax[n,m].imshow(image_bandi, cmap=\"Reds\")\n        ax[n,m].grid(False)\n        band += 1\n        \nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some part of Puerto Rico are not covered by GLDAS.  "},{"metadata":{},"cell_type":"markdown","source":"## Measurements frequency\n\nLet's see what if the overlap in provided data.  \n(inspired from https://www.kaggle.com/jyesawtellrickson/data-measurement-frequency )"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gldas_files = os.listdir(gldas_path)\ngldas_dates = [datetime.strptime(g, 'gldas_%Y%m%d_%H%M.tif') for g in gldas_files]\n\ngfs_files = os.listdir(gfs_path)\ngfs_dates = [datetime.strptime(g, 'gfs_%Y%m%d%H.tif') for g in gfs_files]\n\ns5p_files = os.listdir(sentinel_path)\ns5p_dates = [datetime.strptime(g[:16], 's5p_no2_%Y%m%d') for g in s5p_files]\n\n\nall_dates = (pd.DataFrame(columns=['dataset', 'datetime'])\n             .append(pd.DataFrame(gldas_dates, columns=['datetime'])\n                     .assign(dataset = 'gldas'), sort=True)\n             .append(pd.DataFrame(gfs_dates, columns=['datetime'])\n                     .assign(dataset = 'gfs'), sort=True)\n             .append(pd.DataFrame(s5p_dates, columns=['datetime'])\n                     .assign(dataset = 's5p'), sort=True)\n            ).assign(date = lambda x: x.datetime.apply(lambda x: x.date()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some info on measurements rate:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"all_dates.groupby('dataset').date.agg(\n    min=min,\n    max=max,\n    measurement_period=  lambda x: (x.max()-x.min()).days+1,\n    measurement_count= 'count',\n    measurements_per_day= lambda x: x.count() / ((x.max()-x.min()).days+1)\n).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"daily_data = (pd.date_range('2018-07-01', '2019-06-30').to_frame() # Get the date index to work with\n              .merge( #Merge with all_dates\n                  all_dates.groupby(['dataset', 'date']).date.count().unstack(level=0),\n                  left_index=True,\n                  right_index=True)\n              .drop(columns=[0], axis=1) #remove col used for merge\n             )\n\nsns.heatmap(daily_data.transpose())\nfig = plt.gcf()\nfig.set_size_inches(11,3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"gfs and gldas are available every day. s5p shows missing days."},{"metadata":{},"cell_type":"markdown","source":"# Emission factor estimation\n\nLet's try to see what could be an emission factor of a power plant. First, define one area of 'effect' around each power plant"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# from https://www.kaggle.com/maxlenormand/simplified-emission-for-each-plant-wip\nbuffered_power_plants = global_power_plants_gdf.copy()\nbuffered_power_plants['geometry'] = global_power_plants_gdf.geometry.buffer(0.05)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"lat=18.200178; lon=-66.664513 #puerto rico\nMap = initMap(global_power_plants, lat, lon)   \n\nfor power_plants in range(buffered_power_plants.shape[0]):\n    folium.GeoJson(buffered_power_plants.geometry[power_plants]).add_to(Map)\n\nMap","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To avoid multiple impact, let's choose one power plant away from large cities, and far from others.\n\n![](https://upload.wikimedia.org/wikipedia/en/f/f2/Population_Density%2C_PR%2C_2000_%28sample%29.jpg)\n\nThe oil power plant in north closed to Arecibo might be a good staring point to test some ideas."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"oil_pp_Cambalache = global_power_plants_gdf[global_power_plants_gdf['name'] == 'Cambalache']\noil_pp_Cambalache","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Time series\n\nLet's plot the evolution of N2O over time around our oil power plant.  \nWe'll use the tropospheric_NO2_column_number_density value as it seems to be more related to ground value, and thus emissions from the power plant."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"dates = []\nstats = []\nfor s5p_file in s5p_files:\n    image_path = sentinel_path + s5p_file\n    date = datetime.strptime(s5p_file[:16], 's5p_no2_%Y%m%d')\n    stat = rasterstats.zonal_stats(oil_pp_Cambalache.geometry.to_json(),\n                                   image_path,\n                                   band=2, #2: tropospheric_NO2_column_number_density\n                                   stats=['mean'])\n    stat = stat[0] # get location of pp\n    stat = stat['mean'] # retrieve stat\n    dates.append(date)\n    stats.append(stat)\n\nresults = pd.DataFrame(index=dates, data=stats, columns=['oil'])\nresults.plot()\nplt.title('2: tropospheric_NO2_column_number_density over time in gas_pp_Cambalache')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Difficult to extract tendancies from this plot.  \nLet's use a solar power plant as a comparison. We'll use the one in the north west, as it's close from our test one and far from everything else."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"solar_pp_Oriana = global_power_plants_gdf[global_power_plants_gdf['name'] == 'Oriana Solar Farm']\nsolar_pp_Oriana","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's plot again the n2o for this particular pp."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"dates = []\nstats = []\nfor s5p_file in s5p_files:\n    image_path = sentinel_path + s5p_file\n    date = datetime.strptime(s5p_file[:16], 's5p_no2_%Y%m%d')\n    stat = rasterstats.zonal_stats(solar_pp_Oriana.geometry.to_json(),\n                                   image_path,\n                                   band=2, #2: tropospheric_NO2_column_number_density\n                                   stats=['mean'])\n    stat = stat[0] # get location of pp\n    stat = stat['mean'] # retrieve stat\n    dates.append(date)\n    stats.append(stat)\n\nresults_solar = pd.DataFrame(index=dates, data=stats, columns=['solar'])\nresults_solar.plot()\nplt.title('2: tropospheric_NO2_column_number_density over time in gas_pp_Oriana')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can assume that most of N2O diff between the two pp is due to polution coming from the Oil one"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"results['solar'] = results_solar['solar']\nresults['difference'] = results['oil'] - results['solar']\nresults.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"mean_diff_oil_solar = results['difference'].mean()\nprint(f'Mean difference in tropospheric_NO2_column_number_density over a year: {mean_diff_oil_solar}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**This value might be a starting point to construct an EF based on sp5 ?**  \n\nLet's do the same calculation for all fossil power plants."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fossil_power_plants_gdf = global_power_plants_gdf[(global_power_plants_gdf['primary_fuel'] == 'Oil') |\n                                                  (global_power_plants_gdf['primary_fuel'] == 'Gas') |\n                                                  (global_power_plants_gdf['primary_fuel'] == 'Coal')]\n\nstats_all = []\nfor idx, fossil_power_plant in fossil_power_plants_gdf.iterrows():\n    dates = []\n    stats = []\n    for s5p_file in s5p_files:\n        image_path = sentinel_path + s5p_file\n        date = datetime.strptime(s5p_file[:16], 's5p_no2_%Y%m%d')\n        stat = rasterstats.zonal_stats(fossil_power_plant.geometry,\n                                       image_path,\n                                       band=2, #2: tropospheric_NO2_column_number_density\n                                       stats=['mean'])\n        stat = stat[0] # get location of pp\n        stat = stat['mean'] # retrieve stat\n        dates.append(date)\n        stats.append(stat)\n    name = fossil_power_plant['gppd_idnr']\n    results_fossil_tmp = pd.DataFrame(index=dates, data=stats, columns=[name])\n    stats_all.append(results_fossil_tmp)\n\nresults_fossil = pd.concat(stats_all, axis=1)\nresults_fossil.plot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see the difference compared to the same solar power plant used before."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"all_diff = []\nfor idnr in results_fossil.columns:\n    diff = results_fossil[idnr] - results_solar['solar']\n    #Add more infos from fossil_power_plants_gdf\n    fossil_power_plant = fossil_power_plants_gdf[fossil_power_plants_gdf['gppd_idnr'] == idnr]\n    capacity = float(fossil_power_plant['capacity_mw'])\n    generatedpw = float(fossil_power_plant['estimated_generation_gwh'])\n    fuel = fossil_power_plant['primary_fuel'].item()\n    commissioning_year = int(fossil_power_plant['commissioning_year'])\n    lat = float(fossil_power_plant['lat'])\n    lon = float(fossil_power_plant['lon'])\n    #put infos in a list\n    all_diff.append([idnr, commissioning_year, fuel, capacity, generatedpw, diff.mean(), lat, lon])\n\nall_diff_df = pd.DataFrame(all_diff, columns=['gppd_idnr', 'commissioning_year', 'primary_fuel', 'capacity_mw', 'estimated_generation_gwh', 'EF_N2O', 'lat', 'lon'])\n\ngeometry_power_plants = [Point(x,y) for x,y in zip(all_diff_df['lon'], all_diff_df['lat'])]\nfossil_power_plants_EF_gdf = gpd.GeoDataFrame(all_diff_df, crs = {'init': 'epsg: 4326'}, geometry = geometry_power_plants)\nfossil_power_plants_EF_gdf.head(3)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Plot to see if there is a correlation with fuel type and generated power"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"sns.scatterplot(x=\"estimated_generation_gwh\", y=\"EF_N2O\",\n                hue=\"primary_fuel\", size=\"capacity_mw\",\n                data=all_diff_df[all_diff_df['primary_fuel'] != 'Coal'])\n# control x and y limits\nplt.ylim(-0.000006, 0.000007)\n# Put the legend out of the figure\nplt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Nothing clear.  \nLet's now see how EF evolves with location"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fossil_power_plants_EF_gdf.plot(column='EF_N2O', legend=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have a positive EF only in the north, where urban activities are the highest.  \nLow chance that our EF corresponds to power plant emissions.\n\nMaybe with another point of comparison? Or with taking into account wind ?"},{"metadata":{},"cell_type":"markdown","source":"[WIP]"},{"metadata":{},"cell_type":"markdown","source":"# Google Earth Engine\n\nSome code with GEE to vizualize and access the data in an easier way"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import ee\nfrom kaggle_secrets import UserSecretsClient\nfrom google.oauth2.credentials import Credentials\n\n# Trigger the authentication flow.\n#ee.Authenticate()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#!cat ~/.config/earthengine/credentials","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"user_secret = \"gee\" # 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 GEE\nee.Initialize(credentials=credentials)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Define a method for displaying Earth Engine image tiles to folium map.\ndef add_ee_layer(self, ee_image_object, vis_params, name):\n  map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)\n  folium.raster_layers.TileLayer(\n    tiles = map_id_dict['tile_fetcher'].url_format,\n    attr = \"Map Data © Google Earth Engine\",\n    name = name,\n    overlay = True,\n    control = True\n  ).add_to(self)\n    \n# Add EE drawing method to folium.\nfolium.Map.add_ee_layer = add_ee_layer","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Use this cell to define the collection to see\n\nstartdate = '2018-07-01'\nfinishdate = '2018-07-03'\n\n#Sentinel-5P OFFL NO2: Offline Nitrogen Dioxide\nee_s5p_no2 = (ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')\n              .select('NO2_column_number_density')\n              .filterDate(startdate, finishdate)\n             )\nband_viz_s5p_no2 = {\n    'min': 0,\n    'max': 0.0002,\n    'palette': ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']}\n\n#GFS: Global Forecast System 384-Hour Predicted Atmosphere Data\nee_gfs = (ee.ImageCollection('NOAA/GFS0P25')\n          .select('temperature_2m_above_ground')\n          .filterDate(startdate, finishdate)\n         )\nband_viz_gfs = {\n    'min': -40.0,\n    'max': 35.0,\n    'palette': ['blue', 'purple', 'cyan', 'green', 'yellow', 'red']}\n\n#GLDAS-2.1: Global Land Data Assimilation System\nee_gldas = (ee.ImageCollection('NASA/GLDAS/V021/NOAH/G025/T3H')\n          .select('AvgSurfT_inst')\n          .filterDate(startdate, finishdate)\n         )\nband_viz_gldas = {\n    'min': 250.0,\n    'max': 300.0,\n    'palette': ['1303ff', '42fff6', 'f3ff40', 'ff5d0f']}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Use this cell to add GEE layer to folium\n\n#Init Puerto Rico Map with power plants\nlat=18.200178; lon=-66.664513 #puerto rico\nMap = initMap(global_power_plants, lat, lon)\n\n# Add the sentinel N2O layer to the map object.\nMap.add_ee_layer(ee_s5p_no2.mean(), band_viz_s5p_no2, 'S5P_NO2')\n\n# Add the GFS layer to the map object.\n#Map.add_ee_layer(ee_gfs.mean(), band_viz_gfs, 'GFS')\n\n# Add the GLDAS layer to the map object.\nMap.add_ee_layer(ee_gldas.mean(), band_viz_gldas, 'GLDAS')\n\n# Add a layer control panel to the map.\nMap.add_child(folium.LayerControl())\n\n# Display the map.\nMap","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Next\n\n* Clean code and write functions instead of multiple copy paste (sorry!)\n* Continue to explore EF with time series"}],"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":1}