{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Objective\n\nIn this notebook i will perform a basic exploratory data analysis of the data, also we will make a simple model to calculate the emission factor coefficient of electricity that produce Green House Gases for Puerto Rico.\n\nHere is a link of a public kernel that help me understand better the problem we are trying to solve:\n\nhttps://www.kaggle.com/parulpandey/understanding-the-data-wip\n\n* The model needs to produce a value for the an annual average historical grid-level electricity emissions factor (based on rolling 12-months of data from July 2018 - July 2019) for the sub-national region?\n\nWe also recieve bonuses for other objectives but we will not cover them in this initial notebook"},{"metadata":{},"cell_type":"markdown","source":"# General Equation for Emission Estimation\n\nThe general equation for calculating emission is the following.\n\nE = A X EF X (1 - ER / 100)\n\n* E = emissions\n* A = activity rate\n* EF = emission factor\n* ER = overall emission reduction efficiency %."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport scipy\nimport seaborn as sns\nimport glob\nimport cv2\nfrom tqdm import tqdm_notebook as tqdm\nimport 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\nimport plotly.express as px\npd.set_option('max_columns', 1000)\npd.set_option('max_rows', 1000)\nimport warnings\nwarnings.filterwarnings('ignore')\nimport gc\n\nimport rasterio as rio\nimport folium\nimport tifffile as tiff\n\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn import metrics\nfrom sklearn import preprocessing","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Global Power Plant Database"},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"gpp_df = pd.read_csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\ngpp_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have a lot of usefull data in this dataset, starting by the estimate annual electricity generation in gigawatt-hours.\n\nLet's first check the location of this power plants. Im ussing a function extracted from this great kernel!.\n\nhttps://www.kaggle.com/paultimothymooney/how-to-get-started-with-the-earth-engine-data"},{"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)\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],df[longitude_column].iloc[i]],popup=popup).add_to(plot)\n    return(plot)\n\ndef 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\n\ndef plot_scaled(file_name):\n    vmin, vmax = np.nanpercentile(file_name, (5,95))  # 5-95% stretch\n    img_plt = plt.imshow(file_name, cmap='gray', vmin=vmin, vmax=vmax)\n    plt.show()\n\ndef 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":"gpp_df = split_column_into_new_columns(gpp_df,'.geo','latitude',50,66)\ngpp_df = split_column_into_new_columns(gpp_df,'.geo','longitude',31,48)\ngpp_df['latitude'] = gpp_df['latitude'].astype(float)\na = np.array(gpp_df['latitude'].values.tolist()) # 18 insted of 8\ngpp_df['latitude'] = np.where(a < 10, a + 10, a).tolist()\nlat = 18.200178; lon = -66.664513\nplot_points_on_map(gpp_df, 0, 425, 'latitude', lat, 'longitude', lon, 9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this map we can clearly see that the plant are in Puerto Rico.\n\nLet's check how this plants work, in other words what is the primary fuel for the generation of electricity\n\nAlso we are going to  inspect other variables likes the generation_gwh_2013, generation_gwh_2014 etc..."},{"metadata":{"trusted":true},"cell_type":"code","source":"years = [2013, 2014, 2015, 2016, 2017]\nprint([(gpp_df[f'generation_gwh_{x}'].nunique()) for x in years])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It's look that this columns have only one number (0). Let's continue"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('There are {} power plants'.format(gpp_df.shape[0]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def bar_plot(df, column, title, width, height, n, get_count = True):\n    if get_count == True:\n        cnt_srs = df[column].value_counts()[:n]\n    else:\n        cnt_srs = df\n        \n    trace = go.Bar(\n        x = cnt_srs.index,\n        y = cnt_srs.values,\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 = 'bar_plot')\nbar_plot(gpp_df, 'primary_fuel', 'Primary Fuel Distribution', 800, 500, 100)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Hydro followed by gas are the most common fuels of Puerto Rico electricity plants, coal is the less common\n* We have only one coal plant!"},{"metadata":{"trusted":true},"cell_type":"code","source":"pf_generation = gpp_df.groupby('primary_fuel')['estimated_generation_gwh'].sum()\nbar_plot(pf_generation, 'primary_fuel', 'Electricity Generation Sum by Primary Fuel', 800, 500, 100, False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Coal is the fuel type that produce more energy, on the other hand is the one that power less plants! (1 plant). Coal in known to be an energy commodity\n\nCommissioning year is important, older plants pollute more. Let's check this feature."},{"metadata":{"trusted":true},"cell_type":"code","source":"gpp_df['commissioning_year'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* We have a lot of 0 values, where 0 can be unknown commissioning_year. Also we have a plant that was built in 1942!. \n\nLet'check the capacity mesured in mega watts"},{"metadata":{"trusted":true},"cell_type":"code","source":"pf_capacity = gpp_df.groupby('primary_fuel')['capacity_mw'].sum()\nbar_plot(pf_capacity, 'primary_fuel', 'Capacity Sum by Primary Fuel', 800, 500, 100, False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Oil power plants have the biggest capacity (sum of all oil plants) followed by gas\n\nCapacity and Generation stats"},{"metadata":{"trusted":true},"cell_type":"code","source":"gpp_df.groupby(['primary_fuel']).agg({'estimated_generation_gwh': ['nunique', 'sum', 'mean', 'max', 'min'], 'capacity_mw' : ['nunique', 'sum', 'mean', 'max', 'min']}).reset_index()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Oil fuel power plants have the highest mean capacity, followed by the single Coal power plant"},{"metadata":{},"cell_type":"markdown","source":"Let's check the most important plants (electricity generation)"},{"metadata":{"trusted":true},"cell_type":"code","source":"df = gpp_df[['name','latitude','longitude','primary_fuel','capacity_mw','estimated_generation_gwh', 'owner']].sort_values('estimated_generation_gwh', ascending = False)\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* A.E.S Corp it the plant that generates more energy but's its capacity is not the biggest."},{"metadata":{},"cell_type":"markdown","source":"# Let's Jump to Copernicus image data\n\n* How do we read this data?"},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/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":{},"cell_type":"markdown","source":"Let's check another picture."},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180702T173526_20180708T192358.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":"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":"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":"* Strange, it seems their is no red on this picture"},{"metadata":{"trusted":true},"cell_type":"code","source":"image = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/s5p_no2_20180707T174140_20180713T191854.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":"# Question so Far\n\nPuerto 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.\n\n* NO2 reading are only cause by the power plants electricity generation or more activities? I think they are more factors involved\n* How much pollution of the total pollution does power plants electricity generation produce?\n* How can we isolate the pollution attributable to power plants generation in the remote sensing data?"},{"metadata":{},"cell_type":"markdown","source":"# First Simple Model\n\nLet's calculate the factor with the last image of NO2\n\nWe know that not all energy fuel types produce pollution, the first assumption is going to be that only Caol, Gas and Oil produce pollution.\n\nSearching in google I found that 14% of the pollution is made from power plants electricity generation (im not entirely sure that that is correct).\n\nWe are going to use Simplified emission factor formula:\n\nE / A = EF\n\n* E = emissions\n* A = activity rate\n* EF = emission factor"},{"metadata":{"trusted":true},"cell_type":"code","source":"tiff.imread(image).shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As you can see this photo we have 3 dimensions, the last dimension corresponds to the band. If you look in the data section, their is a description of the bands (last channel)\n\nhttps://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S5P_OFFL_L3_NO2\n\n* NO2 emissions are measured in band1, band2, band3 and band4"},{"metadata":{"trusted":true},"cell_type":"code","source":"p_fuel_types = ['Coal', 'Oil', 'Gas']\n# only consider pollute fuel types\np_type_df = gpp_df[gpp_df['primary_fuel'].isin(p_fuel_types)]\n# sum the electricity generation\np_type_sum = p_type_df['estimated_generation_gwh'].sum()\n# 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\n# use the simplified emission factor formula\nfactor = sum_no2_emission_oe / p_type_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":"# What could be wrong with this simple methology:\n\n* Remember that 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.\n\n* Pictures are different each day, meaning that using only one picture makes no sense at all.\n\n* Consider 14% of pollution is made from power plants electricity? This is a general stat, we should some how use the information given in the power plant database to estimate better. I believe that fuel type and commissioning_year is a very important factor to consider. In other words a really important objective is to calculate the % of pollution made from Puerto Rico power plants!!!.\n\n* Another important factor is to check each plant. Does all plants in Puerto Rico pollute?\n\nLet's try and make better and more accurate baseline factor using all the NO2 emission pictures that we have!!"},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_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":"Let's preprocess this data."},{"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(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":{},"cell_type":"markdown","source":"# NO2 Analysis\n\nFor speed purpose we are going to split the data in 2 dataframes"},{"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":{},"cell_type":"markdown","source":"There most common range between dates is 6 days.\n\nWe have some negatives values, is that correct?\n\nMMMMMM we have some rows were the sum of NO2 emission is NaN. Let's check why"},{"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":"We have some nan values in our images. \n\n* Why do we have some nan values in our images? \n* Should we impute these values for better data quality or ignore them??\n\nLets ignore them for know, going to leave them in my backlog :)"},{"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":"s5p_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":{},"cell_type":"markdown","source":"Now we have a time series dataframe!!. Let's make some plots to visualize the information"},{"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')\nline_plot(s5p_no2_pictures_stats_ig_nan, 'start_date', 'no2_emission_sum', 'NO2 emission by date', 1400, 600)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* We have information of NO2 emission between 2018-07-01 and 2019-06-29\n* We can see a lot of peaks, is there a reasonable explanation for this low values? Maybee ignoring nan values is not correct.\n* Why do we have 387 observations? A year have 365 days\n* We have some duplicate dates"},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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":{"trusted":true},"cell_type":"markdown","source":"* This peaks could be related with missing values!!!\n\n* April 15 of 2019 is an outlier?"},{"metadata":{"trusted":true},"cell_type":"code","source":"line_plot(s5p_no2_pictures_stats[s5p_no2_pictures_stats['start_date']!='2019-04-15'].dropna(), 'start_date', 'no2_emission_sum', 'NO2 emission by date', 1400, 600)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Can we use only the values of the previous graph for our factor?, i think we cant because we are loosing too much information. We need to find a way to deal with this"},{"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":{},"cell_type":"markdown","source":"* Why do we have duplicate days?\n* What should we do with this days?"},{"metadata":{},"cell_type":"markdown","source":"Let's interpolate the nan values!!!"},{"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":{},"cell_type":"markdown","source":"Great! we have clean our NO2 dataset. Let's plot it again."},{"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":{},"cell_type":"markdown","source":"* We have a nice and stable curve\n* Between August 13 of 2018 and March 9 of 2019  we can see a decrease in NO2 emissions (lower trend)\n* What can we do with the missing days? Can we predict them with an arima model?\n\nNow that we have better NO2 information, let's calculate again the factor using the simple methodology"},{"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 / p_type_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":"Our factor changed a lot!!!.\n\nLet's see how we can add more variable to the equation. Let's continue with the weather data!!"},{"metadata":{},"cell_type":"markdown","source":"# Weather Data\n\nWe also have pictures for this data were we can found different bands. Actually their are 9 bands. You can read more in the next link: https://developers.google.com/earth-engine/datasets/catalog/NOAA_GFS0P25"},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/*'\nweather_pictures_path = glob.glob(weather_path)\nlen(weather_pictures_path)\nprint('We have {} pictures of the global forecast system'.format(len(weather_pictures_path)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* We have more than one photo for each day. Probably the last 2 digits before .tif is the hour"},{"metadata":{"trusted":true},"cell_type":"code","source":"tiff.imread(weather_pictures_path[0]).shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Can only see 6 bands in this image. We are going to assume that their are the first 6 bands (plz correct me if i am wrong)"},{"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\nimg_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 img_columns]].isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Weather image dont have missing values, great!"},{"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 img_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":{},"cell_type":"markdown","source":"Let's get the mean to have one observation for each day. Then let's plot them and check the graphs"},{"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 falls between October and June. Let's check the correlation between temperature and NO2 emission!"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Weather data have all the dates, on the other hand some days in the N02 dataframe are missing \nno2_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":"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":"Temperature and relative humidity are the most correlated features with N02"},{"metadata":{},"cell_type":"markdown","source":"# Regression Analysis with Temperature and NO2\n\n* In this experiment we are going make a regression model where we want to predict N02 emission\n\n* We are only going to use weather data, later we can include the electricity produce by the power plants and get the coefficient as our factor!!\n\n* Let's leave the last week as our test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"def preprocess_data(no2_weather, use_lags = True, use_time = True):\n    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']]\n    # get month for groupkfold validation\n    reg_dataset['month'] = reg_dataset['date'].dt.month\n    if use_time:\n        # get day of week as feature\n        reg_dataset['dayofweek'] = reg_dataset['date'].dt.dayofweek\n        # one hot encoder \n        reg_dataset = pd.get_dummies(reg_dataset, columns = ['dayofweek'])\n    if use_lags:\n        # get no2_emissions lags\n        reg_dataset['no2_emission_sum_t1'] = reg_dataset['no2_emission_sum'].shift(1)\n        reg_dataset['no2_emission_sum_t2'] = reg_dataset['no2_emission_sum'].shift(2)\n        reg_dataset['no2_emission_sum_t3'] = reg_dataset['no2_emission_sum'].shift(3)\n        reg_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\n        # drop nan columns produce by the lags\n        reg_dataset.dropna(inplace = True)\n    # split train and test\n    train = reg_dataset[reg_dataset['date'] < '2019-06-23']\n    test = reg_dataset[reg_dataset['date'] >= '2019-06-23']\n    features = [col for col in train.columns if col not in ['date', 'no2_emission_sum', 'month']]\n    return train, test, features\n\n\ndef train_linear_regression(train, test, features, n_folds = 12):\n    # 12 folds, each one representing 1 month\n    target = 'no2_emission_sum'\n    kfold = GroupKFold(n_folds)\n    oof = np.zeros(len(train))\n    predictions = np.zeros(len(test))\n    all_coef = pd.DataFrame()\n    for fold, (trn_ind, val_ind) in enumerate(kfold.split(train[features], groups = train['month'])):\n        print(f'Training and evaluatin fold {fold}')\n        x_train, y_train = train[features].iloc[trn_ind], train[target].iloc[trn_ind]\n        x_val, y_val = train[features].iloc[val_ind], train[target].iloc[val_ind]\n        # standarize train and eval\n        scaler = preprocessing.StandardScaler()\n        x_train = scaler.fit_transform(x_train)\n        x_val = scaler.transform(x_val)\n        test_scaled = scaler.transform(test[features])\n        month = train[['month']].iloc[val_ind]['month'].unique()[0]\n        model = LinearRegression().fit(x_train, y_train)\n        fold_prediction = model.predict(x_val)\n        fold_error = np.sqrt(metrics.mean_squared_error(y_val, fold_prediction))\n        print(f'Our rmse for month {month} is {fold_error}')\n        oof[val_ind] = fold_prediction\n        predictions += model.predict(test_scaled) / n_folds\n        coef = pd.DataFrame({'features': train[features].columns})\n        coef['coef_'] = model.coef_\n        all_coef = pd.concat([all_coef, coef])\n    oof_rmse = np.sqrt(metrics.mean_squared_error(train[target], oof))\n    test_error = np.sqrt(metrics.mean_squared_error(test[target], predictions))\n    fig, ax = plt.subplots(2, 1, figsize = (14, 14))\n    ax[0].plot(train['date'], train[target], color = 'red', label = 'real')\n    ax[0].plot(train['date'], oof, color = 'blue', label = 'prediction')\n    ax[0].set_title('out of fold prediction vs real target')\n    ax[1].plot(test['date'], test[target], color = 'red', label = 'real')\n    ax[1].plot(test['date'], predictions, color = 'blue', label = 'prediction')\n    ax[1].set_title('test prediction vs real target')\n    plt.show()\n    print('The standard deviation for no2 emissions for each month is:')\n    print(train.groupby('month')[target].std().reset_index())\n    print(f'Our out of folds rmse is {oof_rmse}')\n    print(f'Our test rmse is {test_error}')\n    return oof, predictions, all_coef\n\ndef plot_coef(coef):\n    plt.figure(figsize = (12, 8))\n    # absolute for better visuals\n    sns.barplot(abs(coef['coef_']), coef['features'], orient = 'h')\n    plt.title('Feature coefficients')\n    plt.show()\n\n# train with lags and time features\ntrain1, test1, features1 = preprocess_data(no2_weather, use_lags = True, use_time = True)\noof1, predictions1, all_coef1 = train_linear_regression(train1, test1, features1, n_folds = 12)\nplot_coef(all_coef1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"* There are some months that are hardest to predict.\n* Lags are the most important features"},{"metadata":{"trusted":true},"cell_type":"code","source":"train2, test2, features2 = preprocess_data(no2_weather, use_lags = False, use_time = False)\noof2, predictions2, all_coef2 = train_linear_regression(train2, test2, features2, n_folds = 12)\nplot_coef(all_coef2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* If we don't consider lags and time features our model perform badly.\n* Weather data is no that strong to predict NO2 emissions."},{"metadata":{},"cell_type":"markdown","source":"# Still on work!!\n\n* Have a lot of question about the data, need to read more about it so we can understand better the problem.\n\nNext i will see how we can add electriciy generation of the power plants to our lineal regression"}],"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}