{"cells":[{"metadata":{},"cell_type":"markdown","source":"***Objective***\n\nThis Data Science for Good Competition intends to use **remote sensing** techniques to understand Environmental Emissions. Since the whole concept of Satellite Imagery and can be a little overwhelming, this is just an introductory kernel, where I try to explain the various terms and datasets related to satellite Imagery."},{"metadata":{},"cell_type":"markdown","source":"**Problem Statement: Measuring Emissions factors from Satellite Data ?**\n\nWhat is an emission factor?\n\nThe release of GHG into the atmosphere depends mainly on the activity and the product. In order to estimate GHG emissions per unit of available activity, we need to use a factor called emission factor (EF).\n\nFor example: how many  kg of GHG are emitted by 1 kWh of natural gas ?\n\nThus, the emission factor is the sum of emissions of CO2eq of the human activity described as mass unit of CO2eq / reference flows. For example: the EF for the natural gas is the sum of the combustion (0.205 kg CO2eq / kWh ICV) and the upstream (i.e. the production and transport of the gas) (0.0389 kg CO2eq / kWh ICV).\n\n![](https://cdn-images-1.medium.com/max/800/1*3ToZXr2ObHrT5vlTvxU7pg.png) \n\n\nThe general equation for emissions estimation is:\n\n**E = A x EF x (1-ER/100)**\n\nwhere:\n\nE = emissions; A = activity rate; EF = emission factor, and ER =overall emission reduction efficiency, "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":false,"_kg_hide-output":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 seaborn as sns \nimport matplotlib.pyplot as plt \nimport missingno as msno\nimport rasterio as rio\nimport folium\nimport tifffile as tiff \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\nimport os\nfor dirname, _, 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.","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"data = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv'\ndf = pd.read_csv(data)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* get the lat and log from .geo"},{"metadata":{"trusted":true},"cell_type":"code","source":"def split_latnlog_into_new_columns_fromgeo(dataframe,column_to_split,new_column_one,begin_column_one,end_column_one):\n    for i in range(0, len(dataframe)):\n        dataframe.loc[i, new_column_one] = dataframe.loc[i, column_to_split][begin_column_one:end_column_one]\n    return dataframe","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants = split_latnlog_into_new_columns_fromgeo(df,'.geo','latitude',50,66)\npower_plants = split_latnlog_into_new_columns_fromgeo(df,'.geo','longitude',31,48)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **Now lets make some visualization and understand the data :)**"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8,8))\nsns.catplot('primary_fuel', data= power_plants, kind='count', alpha=0.7, height=6, aspect= 3.5)\n\n# Get current axis on current figure\nax = plt.gca()\n\n# Max value to be set\ny_max = power_plants['primary_fuel'].value_counts().max() \n\n# Iterate through the list of axes' patches\nfor p in ax.patches:\n    ax.text(p.get_x() + p.get_width()/5., p.get_height(),'%d' % int(p.get_height()),\n            fontsize=13, color='blue', ha='center', va='bottom')\nplt.title('Frequency plot of primary_fule', fontsize = 20, color = 'black')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **Hydro and gas stands with more than others**\n\n* **Seem like Oli,Gas and the Coal sre polluted Primary_fule**"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8,8))\nsns.catplot('source', data= power_plants, kind='count', alpha=0.7, height=6, aspect= 3.5)\n\n# Get current axis on current figure\nax = plt.gca()\n\n# Max value to be set\ny_max = power_plants['source'].value_counts().max() \n\n# Iterate through the list of axes' patches\nfor p in ax.patches:\n    ax.text(p.get_x() + p.get_width()/5., p.get_height(),'%d' % int(p.get_height()),\n            fontsize=13, color='blue', ha='center', va='bottom')\nplt.title('Frequency plot of source', fontsize = 20, color = 'black')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8,8))\nsns.catplot('owner', data= power_plants, kind='count', alpha=0.7, height=6, aspect= 3.5)\n\n# Get current axis on current figure\nax = plt.gca()\n\n# Max value to be set\ny_max = power_plants['owner'].value_counts().max() \n\n# Iterate through the list of axes' patches\nfor p in ax.patches:\n    ax.text(p.get_x() + p.get_width()/5., p.get_height(),'%d' % int(p.get_height()),\n            fontsize=13, color='blue', ha='center', va='bottom')\nplt.title('Frequency plot of owner ', fontsize = 20, color = 'black')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **lets take a look at the 'estimated_generation_gwh' as it is the estimated annual electricity generation in gigawatt-hours**"},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import stats\n\nsns.distplot(power_plants['estimated_generation_gwh'] , fit=stats.norm);\n\n# Get the fitted parameters used by the function\n(mu, sigma) = stats.norm.fit(power_plants['estimated_generation_gwh'])\nprint( '\\n mu = {:.2f} and sigma = {:.2f}\\n'.format(mu, sigma))\n\n#Now plot the distribution\nplt.legend(['Normal dist. ($\\mu=$ {:.2f} and $\\sigma=$ {:.2f} )'.format(mu, sigma)],\n            loc='best')\nplt.ylabel('Frequency')\nplt.title('estimated_generation_gwh distribution')\n\n#Get also the QQ-plot\nfig = plt.figure()\nres = stats.probplot(power_plants['estimated_generation_gwh'], plot=plt)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **Seem like the distribution is not good lets apply log transformation and look**"},{"metadata":{"trusted":true},"cell_type":"code","source":"#We use the numpy fuction log1p which  applies log(1+x) to all elements of the column\npower_plants[\"estimated_generation_gwh\"] = np.log1p(power_plants[\"estimated_generation_gwh\"])\n\n#Check the new distribution \nsns.distplot(power_plants['estimated_generation_gwh'] , fit=stats.norm);\n\n# Get the fitted parameters used by the function\n(mu, sigma) =stats.norm.fit(power_plants['estimated_generation_gwh'])\nprint( '\\n mu = {:.2f} and sigma = {:.2f}\\n'.format(mu, sigma))\n\n#Now plot the distribution\nplt.legend(['Normal dist. ($\\mu=$ {:.2f} and $\\sigma=$ {:.2f} )'.format(mu, sigma)],\n            loc='best')\nplt.ylabel('Frequency')\nplt.title('estimated_generation_gwh distribution')\n\n#Get also the QQ-plot\nfig = plt.figure()\nres = stats.probplot(power_plants['estimated_generation_gwh'], plot=plt)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **This looks good :)**"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8,8))\nsns.catplot(x=\"primary_fuel\",\n            y=\"capacity_mw\",\n            data=power_plants,\n            jitter=False,\n           )\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **distribution of 'primary_fuel' with the 'capacity_mw'**"},{"metadata":{},"cell_type":"markdown","source":"* **ok now lets chek for the missing values and the columns that should be Droped**"},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_df = df.isnull().sum(axis=0).reset_index()\nmissing_df.columns = ['column_name','missing_count']\nmissing_df = missing_df.ix[missing_df['missing_count']>0]\nmissing_df = missing_df.sort_values(by='missing_count')\n\nind = np.arange(missing_df.shape[0])\nwidth = 0.5\nfig,ax = plt.subplots(figsize=(12,18))\nrects = ax.barh(ind,missing_df.missing_count.values,color='blue')\nax.set_yticks(ind)\nax.set_yticklabels(missing_df.column_name.values, rotation='horizontal')\nax.set_xlabel(\"Count of missing values\")\nax.set_title(\"Number of missing values in each column\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **So Other_fuel1,other_fuel2,other_fuel3 can be droped**"},{"metadata":{"trusted":true},"cell_type":"code","source":"generation_gwh_years = [\"generation_gwh_2013\", \n                        \"generation_gwh_2014\", \n                        \"generation_gwh_2015\", \n                        \"generation_gwh_2016\",\n                        \"generation_gwh_2017\"]\n\npower_plants.loc[:, generation_gwh_years].sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **all this generation_gwh_years can be Droped**"},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_df = power_plants.sort_values('capacity_mw',ascending=False).reset_index()\npower_plants_df=power_plants_df[['name','latitude','longitude','primary_fuel','owner','capacity_mw','estimated_generation_gwh',]]\npower_plants_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **Now as we have Latitude and longitude we check on the map and see what we will get**"},{"metadata":{"trusted":true},"cell_type":"code","source":"world_d= dict(\n   name=list(power_plants['country']),\n    lat=list(power_plants['latitude']),\n   lon=list(power_plants['longitude']),\n   estimated_generation_gwh =list(power_plants['estimated_generation_gwh'])\n)\nworld_data = pd.DataFrame(dict([ (k,pd.Series(v)) for k,v in world_d.items() ]))\nworld_data = world_data.fillna(method='ffill') \n\n\n# create map and display it\nworld_map = folium.Map(location=[18, -66], zoom_start=9)\n\nfor lat, lon, value, name in zip(world_data['lat'], world_data['lon'], world_data['estimated_generation_gwh'], world_data['name']):\n    folium.CircleMarker([lat, lon],\n                        radius=10,\n                        popup = ('<strong>country</strong>: ' + str(name).capitalize() + '<br>'\n                                '<strong>estimated_generation_gwh</strong>: ' + str(value) + '<br>'),\n                        color='red',\n                        \n                        fill_color='red',\n                        fill_opacity=0.7 ).add_to(world_map)\nworld_map","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **Here on the above map if we want to know the  estimated annual electricity generation in gigawatt-hours of the city just clik in the RED cricle :)**"},{"metadata":{},"cell_type":"markdown","source":"* **Now lets move on the other Satellite images data set **\n\nDATASET SELECTION STARTER PACK\n\n* Global Power Plant database by WRI\n* Sentinel 5P OFFL NO2 by EU/ESA/Copernicus\n* Global Forecast System 384-Hour Predicted Atmosphere Data by NOAA/NCEP/EMC\n* Global Land Data Assimilation System by NASA"},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_tif_img_on_map(file_name,lat,lon,zoom):\n    wor_map = folium.Map([lat, lon], zoom_start=zoom)\n    folium.raster_layers.ImageOverlay(\n        image=file_name,\n        bounds = [[18.6,-67.3,],[17.9,-65.2]],\n        colormap=lambda x: (1, 0, 0, x),\n    ).add_to(wor_map)\n    return wor_map","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180708T172237_20180714T190743.tif'\nband = rio.open(image).read(7)\nprint(band.shape)\nvmin, vmax = np.nanpercentile(band, (5,95))  # 5-95% stretch\nimg_plt = plt.imshow(band, cmap='Oranges', vmin=vmin, vmax=vmax)\nplt.show()\nlatitude=18.1429005246921; longitude=-65.4440010699994\nplot_tif_img_on_map(band,latitude,longitude,9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180707T174140_20180713T191854.tif'\nband = rio.open(image).read(7)\nprint(band.shape)\nvmin, vmax = np.nanpercentile(band, (5,95))  # 5-95% stretch\nimg_plt = plt.imshow(band, cmap='Oranges', vmin=vmin, vmax=vmax)\nplt.show()\nlatitude=18.1429005246921; longitude=-65.4440010699994\nplot_tif_img_on_map(band,latitude,longitude,9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/gfs_2018070106.tif'\nband = rio.open(image).read(6)\nvmin, vmax = np.nanpercentile(band, (5,95))  # 5-95% stretch\nimg_plt = plt.imshow(band, cmap='Oranges', vmin=vmin, vmax=vmax)\nplt.show()\nlatitude=18.1429005246921; longitude=-65.4440010699994\nplot_tif_img_on_map(band,latitude,longitude,9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **ok its better to check the property of the TIF images -:)**"},{"metadata":{"trusted":true},"cell_type":"code","source":"from rasterio.plot import show\nimage = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gldas/gldas_20180701_0300.tif'\n\n#load the image\nband = rio.open(image)\nshow(band)\n#All Metadata for the whole raster dataset\nband.meta","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Driver : Data Format\n* dtype : data type\n* width and Height : The dimensions of the image in  **GFS** **image data**  are : 475 X 148\n* count : There are 12 bands in the image\n* crs: Coordinate Reference Systems which refers to the way in which spatial data that represent the earth’s surface.A particular CRS can be referenced by its EPSG code\n* transform : Affine transform (how raster is scaled, rotated, skewed, and/or translated)\n\n***this in gldas image dataset is with 12 channels iamge***"},{"metadata":{},"cell_type":"markdown","source":"the same we qill check for the Global Land Data Assimilation System by NASA(gfs)"},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/gfs_2018070106.tif'\n\n#load the image\nband = rio.open(image)\nshow(band)\n#All Metadata for the whole raster dataset\nband.meta","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* ***this image in gfs image dataset are in 6 channels***"},{"metadata":{},"cell_type":"markdown","source":"the same we qill check for the Sentinel 5P OFFL NO2 by EU/ESA/Copernicus (s5p_no2)"},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180718T173458_20180724T193016.tif'\n\n#load the image\nband = rio.open(image)\nshow(band)\n#All Metadata for the whole raster dataset\nband.meta","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* ***the images in s5p_no2 data set is with 12 channels*** "},{"metadata":{"trusted":true},"cell_type":"code","source":"quantity_of_electricity_generated = power_plants_df['estimated_generation_gwh'][29:30].values\nprint('Quanity of Electricity Generated: ', quantity_of_electricity_generated)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180718T173458_20180724T193016.tif'\n\naverage_no2_emission = [np.average(tiff.imread(image))]\nprint('Average NO2 emissions value: ', average_no2_emission)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"simplified_emissions_factor = float(average_no2_emission/quantity_of_electricity_generated)\nprint('Simplified emissions factor (S.E.F.) for a single power plant on the island of Vieques =  \\n\\n', simplified_emissions_factor, 'S.E.F. units')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"* **this is the just sample Emission Factor cal **\n\n**to find the Emission Factor we have to get all the varablein the formule\n\n**E = A x EF x (1-ER/100)**\n\nwhere:\n\n**E = emissions; A = activity rate; EF = emission factor, and ER =overall emission reduction efficience \n**\n**so to get all this factors i think we need to include more data and see if it’s possible to use remote sensing techniques to better model emissions factors,and to  develop a methodology to calculate an average historical emissions factor for electricity generation in a sub-national region.**\n\n**planning to make a  methodology which is more accurate will soon update**\n\n**will work more and soon will update**\n\n**please upvote if you like that can make me to work more motive :)**"}],"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}