{"cells":[{"metadata":{},"cell_type":"markdown","source":"PROJECT OVERVIEW:\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."},{"metadata":{},"cell_type":"markdown","source":"PROBLEM STATEMENT:\n\nTo 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 Additional bonus points will be awarded for participants that develop methodologies for calculating marginal emissions factors for the sub-national region.\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).\n\nSuch 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).\n\nThe general equation for emissions estimation is:\n\nE = A x EF x (1-ER/100)\n\nwhere:\n\nE = emissions;\n\nA = activity rate;\n\nEF = emission factor, and\n\nER =overall emission reduction efficiency, %"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"## Importing necessary libraries\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\n\n# Analysing datetime\nimport datetime as dt\nfrom datetime import datetime \n\n# Plotting geographical data\nimport folium\nimport rasterio as rio\n\n# File system manangement\nimport os\n\n# Suppress warnings \nimport warnings\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Exploring the dataset**\n\nThe given gppd_120_pr.csv consists of all the power plants which belongs to the Puerto Rico, an unincorporated territory of the United States located in the northeast Caribbean Sea.The latitude of Puerto Rico is 18.200178, and the longitude is -66.664513. The island has been chosen for the analysis since there are fewer confounding factors from nearby areas. Puerto Rico also offers a unique fuel mix and distinctive energy system layout that should make it easier to isolate pollution attributable to power generation in the remote sensing data."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Checking out number of global power plants in Puerto Rico\ngpp_df = pd.read_csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gpp_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Transposing the table of global power plants\ngpp_df.head().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Looking the shape of dataframe\ngpp_df.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas_profiling\npandas_profiling.ProfileReport(gpp_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#to indicate if any value is missing. Any missing values?\ngpp_df.isnull().values.any()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Total missing values for each feature\ngpp_df.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# replacing numerical variables missing values by 0\nglobal_power_plants = gpp_df.fillna(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# List the global power plants which are used primary fuel.\nsns.barplot(x=global_power_plants['primary_fuel'].value_counts().index,y=global_power_plants['primary_fuel'].value_counts())\nplt.title(\"Global power plants which are used Primary fuel\")\nplt.ylabel('Count')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(global_power_plants['commissioning_year'].value_counts(normalize=True))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From this value counts of global power plants, most of the plants did not commission. The power plants are as old as 1915 and the recent belongs to 2012.\n\n\n\nWho owns the power plants? Let's check it out."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.gcf()\nfig.set_size_inches(8,5)\ncolors = [\"aquamarine\", \"plum\", \"orchid\", \"fuchsia\", \"goldenrod\", \"lavender\", \"olive\", \"lime\", \"turquoise\"]\nglobal_power_plants['owner'].value_counts(ascending=True).plot(kind='bar', color = colors)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"https://www.eia.gov/state/analysis.php?sid=RQ#25 Puerto Rico's electricity is supplied by the Puerto Rico Electric Power Authority (PREPA, also known as Autoridad de Energía Eléctrica, or AEE), a government agency that owns the electricity transmission and distribution systems for the main island, Vieques, and Culebra, as well as 80% of the electricity generating capacity\n\nOverall capacity of power plants\n\nWhen it comes to generation capacity, think maximum power output. Capacity is the amount of electricity a generator can produce when it’s running at full blast. This maximum amount of power is typically measured in megawatts (MW) or kilowatts and helps utilities project just how big of an electricity load a generator can handle."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Overall capacity of the plants\noverall_capacity = global_power_plants['capacity_mw'].sum()\nprint('Overall Capacity: '+'{:.4f}'.format(overall_capacity) + ' MW')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"capacity = (global_power_plants.groupby(['primary_fuel'])['capacity_mw'].sum()).to_frame()\ncapacity = capacity.sort_values('capacity_mw',ascending=False)\ncapacity['percentage_of_total'] = (capacity['capacity_mw']/overall_capacity)*100\ncapacity['percentage_of_total'].plot(kind='bar',color=['orange', 'yellow', 'black', 'orange','cyan','blue'])\ncapacity","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the above observation, Oil as primary fuel add up to 68.33% for generation more power in Puerto Rico and 17.97% to Gas. Coal found as fuel to 7.38% and on the other hand, renewable energy generates the power around 6.30%.\n\n**Estimated Power Generation in GigaWatt.**\n\nOn the other hand, estimated power generation refers to the amount of electricity that is produced over period of specific time. That is generally measured in kW/hr, MW/hr or GW/hr."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Overall estimation of the plants\nest_power_gen_gwh = global_power_plants['estimated_generation_gwh'].sum()\nprint('Total Estimated Capacity in Year: '+'{:.4f}'.format(est_power_gen_gwh) + ' MW')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"generation = (global_power_plants.groupby(['primary_fuel'])['estimated_generation_gwh'].sum()).to_frame()\ngeneration = generation.sort_values('estimated_generation_gwh',ascending=False)\ngeneration['percentage_of_total'] = (generation['estimated_generation_gwh']/est_power_gen_gwh)*100\ngeneration['percentage_of_total'].plot(kind='bar',color=['orange', 'yellow', 'black', 'orange','cyan','blue'])\ngeneration","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"More than 97% of Fossil fuel, are mostly used for generation of power in Puerto Rico than renewable energy about 1.03%."},{"metadata":{"trusted":true},"cell_type":"code","source":"def split_column_into_new_columns(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":"def plot_points_on_map(dataframe,begin_index,end_index,latitude_column,latitude_value,longitude_column,longitude_value,zoom):\n    df = dataframe[begin_index:end_index]\n    location = [latitude_value,longitude_value]\n    plot = folium.Map(location=location,zoom_start=zoom, tiles= 'Stamen Terrain')\n    color={ 'Hydro' : 'lightblue', 'Solar' : 'orange', 'Oil' : 'darkblue', 'Coal' : 'black', 'Gas' : 'lightgray', 'Wind' : 'green' }\n    for i in range(0,len(df)):\n        popup = folium.Popup(str(df.primary_fuel[i:i+1]))        \n        folium.Marker([df[latitude_column].iloc[i],\n                       df[longitude_column].iloc[i]],\n                       popup=popup,icon=folium.Icon(icon_color='red',icon ='bolt',prefix='fa',color=color[df.primary_fuel.iloc[i]])).add_to(plot)\n        \n    return(plot)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants = split_column_into_new_columns(global_power_plants,'.geo','latitude',50,66)\nglobal_power_plants = split_column_into_new_columns(global_power_plants,'.geo','longitude',31,48)\nglobal_power_plants['latitude'] = global_power_plants['latitude'].astype(float)\na = np.array(global_power_plants['latitude'].values.tolist()) \nglobal_power_plants['latitude'] = np.where(a < 10, a+10, a).tolist() \n\nlat=18.200178; lon=-66.664513 # Puerto Rico's co-ordinates\nplot_points_on_map(global_power_plants,0,425,'latitude',lat,'longitude',lon,9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants_df = global_power_plants.sort_values('capacity_mw',ascending=False).reset_index()\nglobal_power_plants_df[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now calculating the production capacity of plant on the basis of days, hour and minute."},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants_df = global_power_plants.sort_values('capacity_mw',ascending=False).reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants[\"capacity_factor\"] = global_power_plants[\"estimated_generation_gwh\"]/(global_power_plants[\"capacity_mw\"]*24*365/1000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants[['name','owner','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh', 'capacity_factor']]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the observation, we can see the Capacity Factor is 113 for that coal Power Plant, which is not possible.\n\nAlso, the capacity factor for the hydro, wind and solar Power Plants is 5.2 0.05 & 0.02, which again seems not possible. However, this is much less important since those are very low capacity plants as compared to Oil and Coal Power Plants."},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants['prod_hrs_per_day']=global_power_plants.estimated_generation_gwh*1000/(global_power_plants.capacity_mw*365)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants['prod_hrs_per_hour']=global_power_plants.estimated_generation_gwh*1000/(global_power_plants.capacity_mw*365)*24","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants['prod_hrs_per_minutes']=global_power_plants.estimated_generation_gwh*1000/(global_power_plants.capacity_mw*365)*24*60","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants[['name','owner','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh', 'capacity_factor', 'prod_hrs_per_day', 'prod_hrs_per_hour', 'prod_hrs_per_minutes']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh', 'capacity_factor', 'prod_hrs_per_day','prod_hrs_per_hour','prod_hrs_per_minutes']].groupby(['primary_fuel']).head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants.groupby(['primary_fuel']).agg({'capacity_mw' : ['nunique', 'sum', 'mean', 'max', 'min']}).reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants.groupby(['primary_fuel']).agg({'estimated_generation_gwh': ['nunique', 'sum', 'mean', 'max', 'min']}).reset_index()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh', 'owner']][power_plants_df['owner']=='PREPA'].sort_values('estimated_generation_gwh', ascending = False).groupby(['primary_fuel','capacity_mw']).head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from ast import literal_eval\n\ndef 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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_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)\nglobal_power_plants.drop(columns=['.geo'], inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Checking how many duplicate values are present per columns"},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://stackoverflow.com/questions/9835762/how-do-i-find-the-duplicates-in-a-list-and-create-another-list-with-them\nfeatures = global_power_plants[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh','gppd_idnr','name','owner','primary_fuel','source','url','wepp_id']].columns.values\nunique_max = []\nfor feature in features:\n    values = global_power_plants[feature].value_counts()\n    unique_max.append([feature, values.max(), values.idxmax()])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.transpose((pd.DataFrame(unique_max, columns=['Feature', 'Max duplicates', 'Value'])).sort_values(by = 'Max duplicates', ascending=False).head(50))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants[\"gppd_idnr\"].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"global_power_plants[\"name\"].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"to_drop = [\"generation_gwh_2013\", \"generation_gwh_2014\", \"generation_gwh_2015\", \"generation_gwh_2016\",\"generation_gwh_2017\", \n           \"other_fuel1\",\"other_fuel2\",\"other_fuel3\",\n           \"geolocation_source\",\"year_of_capacity_data\",\n           \"country\", \"country_long\"]\npower_plants_df = global_power_plants.drop(to_drop, axis=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Analysing the images by using the the Rasterio module**\n\nRasterio is a module for reading and writing several different raster formats in Python.\n\nA raster image is an image file format that is defined by a pixel that has one or more numbers associated with it. The number defines the location, size, or color of the pixels. Raster images are commonly .BMP, .GIF, .JPEG, .PNG, and .TIFF files. Today, almost all of the images you see on the Internet and images taken by a digital camera are raster images.\n\n**Exploring the NO2 data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"def overlay_image_on_puerto_rico(file_name,band_layer,lat,lon,zoom):\n    band = rio.open(file_name).read(band_layer)\n    m = folium.Map([lat, lon], zoom_start=zoom)\n    folium.raster_layers.ImageOverlay(\n        image=band,\n        bounds = [[18.6,-67.3,],[17.9,-65.2]],\n        colormap=lambda x: (1, 0, 0, x),\n    ).add_to(m)\n    return m","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import rasterio as rio\nimport folium\nimport tifffile as tiff\n\nimage = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif'\nlatitude=18.1429005246921; longitude=-65.4440010699994\noverlay_image_on_puerto_rico(image,band_layer=7,lat=latitude,lon=longitude,zoom=8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"average_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":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180706T161914_20180712T200737.tif'\noverlay_image_on_puerto_rico(image,band_layer=7,lat=latitude,lon=longitude,zoom=8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From the above observation, there is no red region means there may be nan value."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Likewise you might want to limit the data to only the region of interest\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":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180704T165720_20180710T184641.tif'\noverlay_image_on_puerto_rico(image,band_layer=7,lat=latitude,lon=longitude,zoom=8)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"average_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":"tiff.imread(image).shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"When we look into this shape of the image, that is in three dimension. The last dimension states to the band. https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S5P_OFFL_L3_NO2"},{"metadata":{},"cell_type":"markdown","source":"**Simple factor model**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# only consider pollute fuel types\nfossil_fuels = ['Coal', 'Oil', 'Gas']\nfossil_fuel_df = global_power_plants[global_power_plants['primary_fuel'].isin(fossil_fuels)]\n# sum the electricity generation\nfossil_fuel_sum = fossil_fuel_df['estimated_generation_gwh'].sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# sum the pollution of the last satellite picture\nsum_No2_emission = np.sum(tiff.imread(image)[:, :, 0 : 4])\n# consider 14% of pollution is made from power plants electricity\nsum_No2_emission_oe = sum_No2_emission * 0.14","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# use the simplified emission factor formula\nfactor = sum_No2_emission_oe / fossil_fuel_sum\nprint(f'Simplified emissions factor for Puerto Rico is {factor} mol * h / m^2 * gw')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"the model needs to reproduce an annual average historical grid-level electricity emissions factor (based on rolling 12-months of data from July 2018 - July 2019) for Puerto Rico."},{"metadata":{"trusted":true},"cell_type":"code","source":"import glob\n\nno2_path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/*'\nno2_pictures_path = glob.glob(no2_path)\nlen(no2_pictures_path)\nprint('We have {} pictures of the Copernicus Sentinel'.format(len(no2_pictures_path)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Preprocessing the data"},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm import tqdm_notebook as tqdm\n\n\n# https://www.kaggle.com/ragnar123/exploratory-data-analysis-and-factor-model-idea\n# this function will help us extract the no2 emission data in a tabular way\ndef read_s5p_no2_pictures_data(only_no2_emissions = True):\n    s5p_no2_pictures = []\n    for num, i in tqdm(enumerate(no2_pictures_path), total = 387):\n        temp_s5p_no2_pictures = {'start_date': [], 'end_date': [], 'data': []}\n        temp_s5p_no2_pictures['start_date'] = no2_pictures_path[num][76:84]\n        temp_s5p_no2_pictures['end_date'] = no2_pictures_path[num][92:100]\n        # only no2 emissions\n        if only_no2_emissions:\n            temp_s5p_no2_pictures['data'] = tiff.imread(i)[:, :, 0 : 4]\n            temp_s5p_no2_pictures['no2_emission_sum'] = np.sum(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_mean'] = np.average(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_std'] = np.std(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_max'] = np.max(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_min'] = np.min(tiff.imread(i)[:, :, 0 : 4])\n            s5p_no2_pictures.append(temp_s5p_no2_pictures)\n        # all Copernicus data\n        else:\n            temp_s5p_no2_pictures['data'] = tiff.imread(i)\n            s5p_no2_pictures.append(temp_s5p_no2_pictures)\n    s5p_no2_pictures = pd.DataFrame(s5p_no2_pictures)\n    s5p_no2_pictures['start_date'] = pd.to_datetime(s5p_no2_pictures['start_date'])\n    s5p_no2_pictures['end_date'] = pd.to_datetime(s5p_no2_pictures['end_date'])\n    s5p_no2_pictures.sort_values('start_date', inplace = True)\n    s5p_no2_pictures.reset_index(drop = True, inplace = True)\n    return s5p_no2_pictures\n\ns5p_no2_pictures_df = read_s5p_no2_pictures_data()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"s5p_no2_pictures_stats = s5p_no2_pictures_df[[col for col in s5p_no2_pictures_df.columns if col not in ['data']]]\ns5p_no2_pictures_data = s5p_no2_pictures_df[['data']]\ndel s5p_no2_pictures_df\ns5p_no2_pictures_stats.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def check_arrays(df, row = 1):\n    band1 = pd.DataFrame(s5p_no2_pictures_data['data'][row][:, :, 0])\n    band2 = pd.DataFrame(s5p_no2_pictures_data['data'][row][:, :, 1])\n    band3 = pd.DataFrame(s5p_no2_pictures_data['data'][row][:, :, 2])\n    band4 = pd.DataFrame(s5p_no2_pictures_data['data'][row][:, :, 3])\n    \n    def check_nan(df):\n        df_nan = df.isnull().values.sum()\n        return df_nan\n    \n    band1_nan = check_nan(band1)\n    band2_nan = check_nan(band2)\n    band3_nan = check_nan(band3)\n    band4_nan = check_nan(band4)\n    \n    print('From row {} we have {} nan values for band1'.format(row, band1_nan))\n    print('From row {} we have {} nan values for band2'.format(row, band2_nan))\n    print('From row {} we have {} nan values for band3'.format(row, band3_nan))\n    print('From row {} we have {} nan values for band4'.format(row, band4_nan))\n\n    return band1, band2, band3, band4\n\nband1, band2, band3, band4 = check_arrays(s5p_no2_pictures_data, row = 4)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we are getting nan values in the NO2 emission bands. Ignoring them for just know."},{"metadata":{"trusted":true},"cell_type":"code","source":"# this function ignore nan values from the images\ndef read_s5p_no2_pictures_data_ignore_nan(only_no2_emissions = True):\n    s5p_no2_pictures = []\n    for num, i in tqdm(enumerate(no2_pictures_path), total = 387):\n        temp_s5p_no2_pictures = {'start_date': [], 'end_date': [], 'data': []}\n        temp_s5p_no2_pictures['start_date'] = no2_pictures_path[num][76:84]\n        temp_s5p_no2_pictures['end_date'] = no2_pictures_path[num][92:100]\n        # only no2 emissions\n        if only_no2_emissions:\n            temp_s5p_no2_pictures['data'] = tiff.imread(i)[:, :, 0 : 4]\n            temp_s5p_no2_pictures['no2_emission_sum'] = np.nansum(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_mean'] = np.nanmean(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_std'] = np.nanstd(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_max'] = np.nanmax(tiff.imread(i)[:, :, 0 : 4])\n            temp_s5p_no2_pictures['no2_emission_min'] = np.nanmin(tiff.imread(i)[:, :, 0 : 4])\n            s5p_no2_pictures.append(temp_s5p_no2_pictures)\n        # all Copernicus data\n        else:\n            temp_s5p_no2_pictures['data'] = tiff.imread(i)\n            s5p_no2_pictures.append(temp_s5p_no2_pictures)\n    s5p_no2_pictures = pd.DataFrame(s5p_no2_pictures)\n    s5p_no2_pictures['start_date'] = pd.to_datetime(s5p_no2_pictures['start_date'])\n    s5p_no2_pictures['end_date'] = pd.to_datetime(s5p_no2_pictures['end_date'])\n    s5p_no2_pictures.sort_values('start_date', inplace = True)\n    s5p_no2_pictures.reset_index(drop = True, inplace = True)\n    return s5p_no2_pictures\n\ns5p_no2_pictures_df_ig_nan = read_s5p_no2_pictures_data_ignore_nan()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://stackoverflow.com/questions/14940743/selecting-excluding-sets-of-columns-in-pandas\ns5p_no2_pictures_stats_ig_nan = s5p_no2_pictures_df_ig_nan[[col for col in s5p_no2_pictures_df_ig_nan.columns if col not in ['data']]]\ndel s5p_no2_pictures_df_ig_nan\ns5p_no2_pictures_stats_ig_nan.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Using plotly.express\nimport plotly.express as px\nfig = px.line(s5p_no2_pictures_stats_ig_nan, x='start_date', y='no2_emission_sum')\nfig.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfrom plotly import tools, subplots\nimport plotly.offline as py\npy.init_notebook_mode(connected = True)\nimport plotly.graph_objs as go\n\ndef line_plot_check_nan(df1, df2, x, y, title, width, height):\n    \n    trace1 = go.Scatter(\n        x = df1[x],\n        y = df1[y],\n        mode='lines',\n        name='with_nans',\n        marker = dict(\n            color = '#1E90FF', \n        ), \n    )\n    \n    df3 = df2.dropna()\n    trace2 = go.Scatter(\n        x = df3[x],\n        y = df3[y],\n        mode='markers',\n        name='no_nans',\n        marker = dict(\n            color = 'red', \n        ), \n    )\n    \n    layout = go.Layout(\n        title = go.layout.Title(\n            text = title,\n            x = 0.5\n        ),\n        font = dict(size = 14),\n        width = width,\n        height = height,\n    )\n    \n    data = [trace1, trace2]\n    fig = go.Figure(data = data, layout = layout)\n    py.iplot(fig, filename = 'line_plot')\nline_plot_check_nan(s5p_no2_pictures_stats_ig_nan, s5p_no2_pictures_stats, 'start_date', 'no2_emission_sum', 'NO2 emission by date', 1400, 600)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From this graph, this peak might be the missing values in images. 15 April, 2019 is an outlier among in the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"n_duplicates_dates = s5p_no2_pictures_stats_ig_nan.shape[0] - s5p_no2_pictures_stats_ig_nan.drop_duplicates(subset = ['start_date', 'end_date']).shape[0]\nprint(f'We have {n_duplicates_dates} duplicate days')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# this function will help us extract the no2 emission data in a tabular way\ndef read_s5p_no2_pictures_data_fill(only_no2_emissions = True):\n    s5p_no2_pictures = []\n    for num, i in tqdm(enumerate(no2_pictures_path), total = 387):\n        temp_s5p_no2_pictures = {'start_date': [], 'end_date': [], 'data': []}\n        temp_s5p_no2_pictures['start_date'] = no2_pictures_path[num][76:84]\n        temp_s5p_no2_pictures['end_date'] = no2_pictures_path[num][92:100]\n        # only no2 emissions\n        if only_no2_emissions:\n            image = tiff.imread(i)[:, :, 0 : 4]\n            band1 = pd.DataFrame(image[: ,: , 0]).interpolate()\n            band1.fillna(band1.mean(), inplace = True)\n            band2 = pd.DataFrame(image[: ,: , 1]).interpolate()\n            band2.fillna(band2.mean(), inplace = True)\n            band3 = pd.DataFrame(image[: ,: , 2]).interpolate()\n            band3.fillna(band3.mean(), inplace = True)\n            band4 = pd.DataFrame(image[: ,: , 3]).interpolate()\n            band4.fillna(band4.mean(), inplace = True)\n            image = np.dstack((band1, band2, band3, band4))\n            temp_s5p_no2_pictures['data'] = image\n            temp_s5p_no2_pictures['no2_emission_sum'] = np.sum(image)\n            temp_s5p_no2_pictures['no2_emission_mean'] = np.average(image)\n            temp_s5p_no2_pictures['no2_emission_std'] = np.std(image)\n            temp_s5p_no2_pictures['no2_emission_max'] = np.max(image)\n            temp_s5p_no2_pictures['no2_emission_min'] = np.min(image)\n            s5p_no2_pictures.append(temp_s5p_no2_pictures)\n        # all Copernicus data\n        else:\n            temp_s5p_no2_pictures['data'] = tiff.imread(i)\n            s5p_no2_pictures.append(temp_s5p_no2_pictures)\n    s5p_no2_pictures = pd.DataFrame(s5p_no2_pictures)\n    s5p_no2_pictures['start_date'] = pd.to_datetime(s5p_no2_pictures['start_date'])\n    s5p_no2_pictures['end_date'] = pd.to_datetime(s5p_no2_pictures['end_date'])\n    s5p_no2_pictures.sort_values('start_date', inplace = True)\n    s5p_no2_pictures.reset_index(drop = True, inplace = True)\n    return s5p_no2_pictures\n\ns5p_no2_pictures_df_fill = read_s5p_no2_pictures_data_fill()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"s5p_no2_pictures_stats_fill = s5p_no2_pictures_df_fill[[col for col in s5p_no2_pictures_df_fill.columns if col not in ['data']]]\ndel s5p_no2_pictures_df_fill\ns5p_no2_pictures_stats_fill.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# drop nan values and check again for duplicate columns\ns5p_no2_pictures_stats_fill = s5p_no2_pictures_stats_fill[s5p_no2_pictures_stats_fill['start_date']!='2019-04-15'].dropna()\n# drop 2019-04-15 (probably an outlier or a rare event that can affect our factor calculation)\nduplicate_columns = s5p_no2_pictures_stats_fill.shape[0] - s5p_no2_pictures_stats_fill.drop_duplicates(subset = ['start_date', 'end_date']).shape[0]\nprint(f'We have {duplicate_columns} duplicate columns')\nprint('We have {} days of data'.format(s5p_no2_pictures_stats_fill['start_date'].nunique()))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def line_plot(df, x, y, title, width, height):\n    trace = go.Scatter(\n        x = df[x],\n        y = df[y],\n        mode='lines',\n        name='lines',\n        marker = dict(\n            color = '#1E90FF', \n        ), \n    )\n    \n    layout = go.Layout(\n        title = go.layout.Title(\n            text = title,\n            x = 0.5\n        ),\n        font = dict(size = 14),\n        width = width,\n        height = height,\n    )\n    \n    data = [trace]\n    fig = go.Figure(data = data, layout = layout)\n    py.iplot(fig, filename = 'line_plot')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"line_plot(s5p_no2_pictures_stats_fill, 'start_date', 'no2_emission_sum', 'NO2 emission by date', 1400, 800)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get the mean NO2 emission between 2018/07/01 and 2019/06/29\nsum_no2_emission = s5p_no2_pictures_stats_fill['no2_emission_sum'].mean()\n# consider 14% of pollution is made from power plants electricity\nsum_no2_emission_oe = sum_no2_emission * 0.14\n# use the simplified emission factor formula (sum of estimated generation from Caol, Oil and Gas plants)\nfactor = sum_no2_emission_oe / fossil_fuel_sum\nprint(f'Simplified emissions factor for Puerto Rico is {factor} mol * h / m^2 * gw')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Emission factor slightly dropped down the previous one if we treat on nan and duplicate values which found in the weather data."},{"metadata":{},"cell_type":"markdown","source":"**Exploring the weather data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_path_data = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/*'\nweather_pictures_path = glob.glob(weather_path_data)\nlen(weather_pictures_path)\nprint('We have {} pictures of the global forecast system'.format(len(weather_pictures_path)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tiff.imread(weather_pictures_path[0]).shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"from the above dimension of images of weather, there are 6 bands.\n\n**Preprocess the weather data:**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# this function will help us extract weather pictures in a tabular way\ndef read_weather_data():\n    weather_pictures = []\n    for num, i in tqdm(enumerate(weather_pictures_path), total = len(weather_pictures_path)):\n        temp_weather_pictures = {'date': [], 'temperature_2m_above_ground': [], 'specific_humidity_2m_above_ground': [], 'relative_humidity_2m_above_ground': [], \n                                 'u_component_of_wind_10m_above_ground': [], 'v_component_of_wind_10m_above_ground': [], 'total_precipitation_surface': []}\n        temp_weather_pictures['date'] = weather_pictures_path[num][68:-6]\n        temp_weather_pictures['date'] = weather_pictures_path[num][68:-6]\n        image = tiff.imread(i)\n        temp_weather_pictures['temperature_2m_above_ground'] = image[ : , : , 0]\n        temp_weather_pictures['specific_humidity_2m_above_ground'] = image[ : , : , 1]\n        temp_weather_pictures['relative_humidity_2m_above_ground'] = image[ : , : , 2]\n        temp_weather_pictures['u_component_of_wind_10m_above_ground'] = image[ : , : , 3]\n        temp_weather_pictures['v_component_of_wind_10m_above_ground'] = image[ : , : , 4]\n        temp_weather_pictures['total_precipitation_surface'] = image[ : , : , 5]\n        temp_weather_pictures['temperature_2m_above_ground_mean'] = np.average(image[ : , : , 0])\n        temp_weather_pictures['specific_humidity_2m_above_ground_mean'] = np.average(image[ : , : , 1])\n        temp_weather_pictures['relative_humidity_2m_above_ground_mean'] = np.average(image[ : , : , 2])\n        temp_weather_pictures['u_component_of_wind_10m_above_ground_mean'] = np.average(image[ : , : , 3])\n        temp_weather_pictures['v_component_of_wind_10m_above_ground_mean'] = np.average(image[ : , : , 4])\n        temp_weather_pictures['total_precipitation_surface_mean'] = np.average(image[ : , : , 5])\n        \n        weather_pictures.append(temp_weather_pictures)\n    \n    weather_pictures = pd.DataFrame(weather_pictures)\n    weather_pictures['date'] = pd.to_datetime(weather_pictures['date'], infer_datetime_format  = True)\n    weather_pictures.sort_values('date', inplace = True)\n    weather_pictures.reset_index(drop = True, inplace = True)\n    return weather_pictures\n\nweather_pictures_df = read_weather_data()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_pictures_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# check missing values\nimage_columns = ['temperature_2m_above_ground', 'specific_humidity_2m_above_ground', 'relative_humidity_2m_above_ground', \n               'u_component_of_wind_10m_above_ground', 'v_component_of_wind_10m_above_ground', 'total_precipitation_surface']\nweather_pictures_df[[col for col in weather_pictures_df.columns if col not in image_columns]].isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Weather image data do not have missing values."},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_pictures_df_stats = weather_pictures_df[[col for col in weather_pictures_df.columns if col not in image_columns]]\nn_duplicates = weather_pictures_df_stats.shape[0] - weather_pictures_df_stats['date'].nunique()\nprint(f'We have {n_duplicates} observations that belongs to a date with one or more records')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_pictures_df_stats = weather_pictures_df_stats.groupby('date').mean().reset_index()\nprint('We have data for {} days'.format(weather_pictures_df_stats['date'].nunique()))\nprint('Our data start on {} and finish in {}'.format(weather_pictures_df_stats['date'].min(), weather_pictures_df_stats['date'].max()))\nline_plot(weather_pictures_df_stats, 'date', 'temperature_2m_above_ground_mean', 'Temperature by Date', 1400, 800)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Temperature is falling down in the month of October 2018. Need to check the correlation in between temperature and No2."},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_weather = s5p_no2_pictures_stats_fill[['start_date', 'no2_emission_sum']].merge(weather_pictures_df_stats, left_on = 'start_date', right_on = 'date', how = 'left')\nno2_tem_corr = no2_weather[['no2_emission_sum', 'temperature_2m_above_ground_mean']].corr().loc['no2_emission_sum', 'temperature_2m_above_ground_mean']\nprint(f'NO2 and temeprature have a correlation of: {no2_tem_corr}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_weather.columns","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"line_plot(weather_pictures_df_stats, 'date', 'specific_humidity_2m_above_ground_mean', 'Specific Humidity by Date', 1400, 800)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"line_plot(weather_pictures_df_stats, 'date', 'relative_humidity_2m_above_ground_mean', 'Relative Humidity by Date', 1400, 800)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"line_plot(weather_pictures_df_stats, 'date', 'u_component_of_wind_10m_above_ground_mean', 'U Component of Wind by Date', 1400, 800)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"line_plot(weather_pictures_df_stats, 'date', 'v_component_of_wind_10m_above_ground_mean', 'V Component of Wind by Date', 1400, 800)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"line_plot(weather_pictures_df_stats, 'date', 'total_precipitation_surface_mean', 'Total Precipitation Surface by Date', 1400, 800)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (14, 8))\nsns.heatmap(no2_weather.corr(), annot = True, cmap = 'coolwarm')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From heatmap, temerature and No2 are correlated with each other."},{"metadata":{"trusted":true},"cell_type":"code","source":"reg_dataset = no2_weather[['date', 'temperature_2m_above_ground_mean', 'specific_humidity_2m_above_ground_mean', 'relative_humidity_2m_above_ground_mean', 'u_component_of_wind_10m_above_ground_mean', \n                               'v_component_of_wind_10m_above_ground_mean', 'total_precipitation_surface_mean', 'no2_emission_sum']]\nreg_dataset['month'] = reg_dataset['date'].dt.month\nreg_dataset['dayofweek'] = reg_dataset['date'].dt.dayofweek\n\nreg_dataset = pd.get_dummies(reg_dataset, columns = ['dayofweek'])\nreg_dataset['no2_emission_sum_t1'] = reg_dataset['no2_emission_sum'].shift(1)\nreg_dataset['no2_emission_sum_t2'] = reg_dataset['no2_emission_sum'].shift(2)\nreg_dataset['no2_emission_sum_t3'] = reg_dataset['no2_emission_sum'].shift(3)\nreg_dataset['no2_emission_rolling_mean_t1t3'] = (reg_dataset['no2_emission_sum_t1'] + reg_dataset['no2_emission_sum_t2'] + reg_dataset['no2_emission_sum_t3']) / 3","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reg_dataset.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reg_dataset.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reg_dataset.fillna(reg_dataset.mean()).head(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install dabl","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import dabl\nfrom dabl import plot\n\nplot(reg_dataset, 'no2_emission_sum')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ec = dabl.SimpleClassifier(random_state=0).fit(reg_dataset, target_col=\"month\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reg_dataset_df = reg_dataset[['date','no2_emission_sum']]\n\n#Set Date column as the index column.\nreg_dataset_df.set_index('date', inplace=True)\nreg_dataset_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15, 5))\nplt.ylabel(\"no2_emission_sum\")\nplt.xlabel(\"date\")\nplt.plot(reg_dataset_df)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Please upvote this notebook. If any suggestion, let me know. Thank you....."}],"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}