{"cells":[{"metadata":{},"cell_type":"markdown","source":"# A simplified emissions factor for each plant - Work in Progress\n\nThis notebook is a very earlier work on this analysis competition, and I hope will grow over time. I also would like to use this notebook as an opportunity to showcase some of the remote sensing tools and methods that I have come across over the months that I've been doing this. I still have a lot to learn, but also know that I'm probably more familiar with remote sensing than most people here on Kaggle.\n\nUpdate 1 (notebook version 4 + 5): Created Sentinel 5 time-series dataframe for a buffered area around each power plant."},{"metadata":{},"cell_type":"markdown","source":"First of all, we start by importing a library that is going to become really handy later on to get some information only around the plants"},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install rasterstats --quiet","execution_count":null,"outputs":[]},{"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 os\nfrom tqdm import tqdm\n\n# Plotting libraries\nimport matplotlib.pyplot as plt\nimport seaborn as sns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Geospatial libraries that we will be using for this\nimport rasterio\nimport rasterstats\nimport folium\nimport geopandas as gpd\nfrom shapely.geometry import Point","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Creating a buffer around the power plants"},{"metadata":{},"cell_type":"markdown","source":"I'll be using the power plants converted into a GeoJSON (accessible in [this dataset](https://www.kaggle.com/maxlenormand/geolocated-power-plants-geojson) from [this notebook](https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson))."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"power_plants_path = '/kaggle/input/geolocated-power-plants-geojson/Geolocated_gppd_120_pr.geojson'\npower_plants = gpd.read_file(power_plants_path)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are only interested in the area around each plant. In order to have only local information for each plant, we can buffer them.\n\nNote here that the unit is in arc degreees, not meters!"},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"buffered_power_plants = power_plants\nbuffered_power_plants['geometry'] = power_plants.geometry.buffer(0.1)","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Plot on the map\nlat=18.200178; lon=-66.664513\n\ndef plot_polygons_on_Map(geodataframe,\n                       original_latitude = lat,\n                       original_longitude = lon,\n                       zoom=9):\n    plot = folium.Map(location = (original_latitude, original_longitude), zoom_start=zoom)\n\n\n    geojson = folium.GeoJson(geodataframe,\n                            style_function=lambda x: {'Color':'red',\n                                                      'fillColor':'yellow' if x['properties']['primary_fuel']=='Solar'  else 'red'})\n    popup = folium.Popup(str(geodataframe.primary_fuel))\n    popup.add_to(geojson)\n    geojson.add_to(plot)\n\n    return(plot)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Yellow filled cirlces are solar plants, the rest are red (Simply to try out some folium libraries so far)."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot_polygons_on_Map(buffered_power_plants)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Looking into all the Sentinel 5 images at our disposal"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Make a dataframe containing all the S5 raster paths and dates\n\nfrom datetime import datetime\n\ndf_s5 = pd.DataFrame()\n\nfiles=[]\ncaptured_datetime = []\nfor dirname, _, filenames in os.walk('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2'):\n    for filename in filenames:\n        files.append(os.path.join(dirname, filename))\n        captured_datetime.append(datetime.strptime(filename.split('_')[-2],'%Y%m%dT%H%M%S').date())\n        \ndf_s5['S5_path'] = files\ndf_s5['Captured_datetime'] = captured_datetime\n\n# Sort dataframe by ascending time\ndf_s5.sort_values('Captured_datetime', inplace=True)\ndf_s5.reset_index(inplace=True, drop=True)\n\ndf_s5.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f\"There are {df_s5.shape[0]} Sentinel 5 images from {df_s5['Captured_datetime'][0]} to {df_s5['Captured_datetime'][df_s5.shape[0]-1]}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Zonal N02 statistics around each power plant\n\nNext step is to get the zonal statistics of the NO2 for each polygon. This is where [rasterstats' zonal statistics come into play](https://pythonhosted.org/rasterstats/manual.html#raster-data-sources)"},{"metadata":{"trusted":true},"cell_type":"code","source":"test_S5_raster_path = df_s5['S5_path'][0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can only use 1 band at a time for the zonal statistics, so we need to know which one to use"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def return_bands(file_name):\n    # Function from the organizers: https://www.kaggle.com/paultimothymooney/explore-image-metadata-s5p-gfs-gldas\n    image = rasterio.open(file_name)\n    for i in image.indexes:\n        desc = image.descriptions[i-1]\n        print(f'{i}: {desc}')\n        \nreturn_bands(test_S5_raster_path)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Just as this is a test to develop the method more than anything, I'll use the 1st band, the N02 column number density. This seems to be what is used in the code provided by the organizers, for now this will be good enough"},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_stats = rasterstats.zonal_stats(power_plants_path,\n                                                 test_S5_raster_path,\n                                                 band=1,\n                                                 stats=['mean'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see what the first few mean N02 around each plant look like:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Plotting the first few elements\npower_plants_stats[:5]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Simplified emissions factor\n\nWe're going to start by saying our simplified emissions factor is: SEF = [N02 measurement] / [Quantity of electricity generated]****"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"cols_to_keep = ['system:index', 'capacity_mw', 'commissioning_year', 'estimated_generation_gwh', 'primary_fuel',  'geometry']\nsimplified_emissions_factor_df = buffered_power_plants.copy(deep=True)\nsimplified_emissions_factor_df = simplified_emissions_factor_df[cols_to_keep]\n\nsimplified_emissions_factor_df['mean_N02'] = [x['mean'] for x in power_plants_stats]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"simplified_emissions_factor_df['Simpl_Emiss_Factor'] = simplified_emissions_factor_df['mean_N02'] / simplified_emissions_factor_df['estimated_generation_gwh']","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"MULTIPLICATION_FACTOR = 100000\n\nplt.figure()\nsns.distplot(simplified_emissions_factor_df['Simpl_Emiss_Factor']*MULTIPLICATION_FACTOR,\n             kde=False, \n             label=f'SEF x {MULTIPLICATION_FACTOR}')\nplt.legend()\nplt.title('Simplified Emission Factor for each plant (N02 / Estimated Emsission)')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As many have pointed out, there seems to be a big outlier. I knew this from the get-go, but just wanted to have something working roughly before tackling that problem."},{"metadata":{},"cell_type":"markdown","source":"## Getting a N02 time-series for each power plant\n\nNow that we have a time-series of Sentinel 5, and a buffer around each power plant, we can get the N02 readings for each Sentinel 5 image, for each power plant\n\nAs a first step, we will only do this over the 1st band, once again."},{"metadata":{},"cell_type":"markdown","source":"This next step is going to take a while to compute (~12 minutes). It could probably be optimized, but that will be for later"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"SUBSET_TO_TEST = df_s5.shape[0]\n\nN02_measurments_df = pd.DataFrame({'index_power_plant':power_plants['system:index']})\n\nfor s5_image in tqdm(range(SUBSET_TO_TEST)):\n    S5_raster_stats = rasterstats.zonal_stats(power_plants_path,\n                                                     df_s5['S5_path'][s5_image],\n                                                     band=1,\n                                                     stats=['mean'])\n    list_of_S5_raster_stats = [plant['mean'] for plant in S5_raster_stats]\n    \n    N02_measurments_df[df_s5['Captured_datetime'][s5_image]] = list_of_S5_raster_stats","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"N02_measurments_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Investigating time-series of 1 power plant\n\nFirst of all, we're just going to look at the N02 measurement of a single plant"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(f'There are {N02_measurments_df.iloc[0].isna().sum()} NaNs ({np.round((N02_measurments_df.iloc[0].isna().sum() / N02_measurments_df.shape[1]) * 100, 2)}% of data)')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Removing NaNs\n\ntest_power_plant_time_series = N02_measurments_df.iloc[0].dropna()\ntest_power_plant_time_series.isna().sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(15, 10))\nplt.plot(list(test_power_plant_time_series[1:]))\nplt.title(f'N02 measurement of plant {N02_measurments_df[\"index_power_plant\"][0]}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Looking at differences between different kinds of power plants\n\nWe'll look if there are any differences for N02 emissions for different types of power plants"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Making a dataframe only keeping the info related to the enegery production\n\npower_plant_type_infos = power_plants[['capacity_mw', 'commissioning_year', 'estimated_generation_gwh', 'primary_fuel']]\npower_plant_type_infos.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Different types of power plants and their counts:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"np.unique(power_plant_type_infos.primary_fuel, return_counts=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"id_coal = power_plant_type_infos.loc[power_plant_type_infos['primary_fuel']=='Coal'].index.values[0]\nid_coal","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# Looking for coal and solar\n\n# --- Coal ----\n\n# Removing NaNs\nid_coal = power_plant_type_infos.loc[power_plant_type_infos['primary_fuel']=='Coal'].index.values[0]\ncoal_plants_stats = N02_measurments_df.iloc[id_coal].dropna()\nprint(f'There are {coal_plants_stats.isna().sum()} NaNs left')\n\n# --- Solar ----\n\n# Removing NaNs (only keeping 1 of the solar plants)\nid_solar = power_plant_type_infos.loc[power_plant_type_infos['primary_fuel']=='Solar'].index.values[0]\nsolar_plants_stats = N02_measurments_df.iloc[id_solar].dropna()\nprint(f'There are {solar_plants_stats.isna().sum()} NaNs left')\n\n\nplt.figure(figsize=(30, 10))\n\nplt.plot(list(coal_plants_stats[1:]), label='Coal', color='brown')\nplt.plot(list(solar_plants_stats[1:]), label='Solar', color='blue')\n\nplt.legend()\nplt.title(f'N02 measurement of the only Coal and 1 solar plant')\nplt.show()\n\n\nplants_of_interest = buffered_power_plants.iloc[[id_solar, id_coal]]\nplot_polygons_on_Map(geodataframe=plants_of_interest)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Conclusions (for now) + future work\n\nThis is a hyper simple approach, and just a small improvement on what the organizers had created. I just wanted to show how rasterstats could be used in this context to get zonal statistics around each plant.\n\nThere are a lot of issues at hand:\n- Here I only used 1 image, while we would want to use all the images of a given year.\n- I also only used the estimated emissions column, but we have no info on when that is, nor for how long? This is where I think things start getting really tricky.\n- For now, there is basically no model, it's just some computing of a very simplified emissions factor. Ideally, we want to be able to know our ground truth, which also seems like it's going to be tricky.\n- We need to find how to get the activation rate, or even the effficiency of the plants somehow to get a more accurate model.\n- As pointed out in some other topcis / notebooks, Sentinel 5 meaures the N02 over the country, wherever it might come from. So we're going to have to find a way to remove all external sources (like cars, housing, etc.) even if it seems like some of the plants might be further from cities than others.\n- There's probably some pre-processing to be done on some of the images. I've never worked on Sentinel 5, but want to look deaper into it.\n- Find more info about weather data, impact of temperature, winds, etc. on N02. (maybe when it's colder, there is more demande for electricity (heating) which would be seen in the N02 emissions?).\n- A lot of other things!"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}