{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":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\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 numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport warnings\nwarnings.filterwarnings('ignore')\nimport os\nimport folium\nimport rasterio as rio\nimport tifffile as tiff\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn import metrics\nfrom sklearn.metrics import mean_squared_error\n\nprint(os.listdir('/kaggle/input/ds4g-environmental-insights-explorer/eie_data'))\n# Any results you write to the current directory are saved as output.\neie_data_path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Global Power Plant Database\n\nThe island of Puerto Rico has multiple power plants for its energy generation. We will begin by seeing where these power plants are."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"gpp_df = pd.read_csv(eie_data_path+'/gppd/gppd_120_pr.csv')\ngpp_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#code source: https://www.kaggle.com/paultimothymooney/overview-of-the-eie-analytics-challenge\ndef 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\n\ngpp_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":{"trusted":true},"cell_type":"code","source":"print('There are ', gpp_df.shape[0],' power plants')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gpp_df.head()\n#For the purposes of the model, we will use estimated power generation in gwh as a feature\ntotal_power_generation_2017 = gpp_df['estimated_generation_gwh'].sum()\nprint('The total estimated power generation for all power plants in puerto for 2017 rico is: ', total_power_generation_2017, ' gwh')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sentinel-5 Precursor Data\n\nThe Sentinel-5 Precuror satellite has sensors to gauge NO2 concentrations. Below we will visualize the NO2 concentrations and display them over the island of Puerto Rico. "},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_emissions_image = eie_data_path+'/s5p_no2/s5p_no2_20180701T161259_20180707T175356.tif'\nlatitude=18.1429005246921; longitude=-65.4440010699994\noverlay_image_on_puerto_rico(no2_emissions_image,band_layer=1,lat=latitude,lon=longitude,zoom=8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the image above, we are overlaying the first band of an image taken on July 7th, 2017 over the island of Puero Rico. The overlayed image (in red) displays the total vertical column of NO 2 in mol/m^2. The shade of red varies based on the quantity. Below there is a similar overlay for band 2 (measurements of the vertical column of NO 2 for the troposphere)."},{"metadata":{"trusted":true},"cell_type":"code","source":"overlay_image_on_puerto_rico(no2_emissions_image,band_layer=2,lat=latitude,lon=longitude,zoom=8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### First Simple Model\n\nAn Emissions Factor is a value that relates the quanity of a pollutant in the atmosphere with an activity associated with the release of that pollutant. \n\nWe can calculate emissions using the following equation: \n\nE = A * EF * (1 - ER / 100 )\n\nHere A refers to an activity, EF is the emissions factor, and ER is the overall emissions reduction efficiency (percentage value). \n\nIf we solve for the Emissions Factor:\n\nEF = (E * 100) / (A * (1-ER))\n\nBy solving for the EF, we can try to calculate the factor for the power plant data\n\nFrom the Sentinel Satellite, NO2 emissions are measured in band1, band2, band3 and band4. The total emissions (E) for a given day (each tif file for emissions is a snapshot for a single day), can be calculated by summing these four bands).\n\nThe Activity for which we're interested in is the energy generation from power plants. One way to calculate the emissions factor for the energy generation would be to use the activity from all the power plants as A, and assume that the emissions reduction efficiency is 0. Therefore:\n\nEF = (E) / (A)"},{"metadata":{"trusted":true},"cell_type":"code","source":"activity = total_power_generation_2017\n#Using the emissions snapshot from above\nno2_emissions_sum = tiff.imread(no2_emissions_image)[:,:,0:4].sum()\nemissions_factor = (no2_emissions_sum ) / (total_power_generation_2017)\nprint(\"Emissions Factor for power plant energy generation activity : \", emissions_factor)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have calculated the emissions factor to be 3.886 * 10^-5. This is a very rough estimate of the emissions factor. This model makes many assumptions:\n\n- It does not take any spatial information into account, i.e. the original NO2 emissions data has measurements over the entire island of Puerto Rico. The emissions (as seen in the overlayed map), vary based on the location. Instead emissions over the entire island are summed. Similar to the NO2 emissions, the power plants total power generation does not take spatial information into account (the activity is not uniform over the entire island). Due to these assumptions, it is assumed that only energy generation activity from power plants contributes to the generation of NO2 emissions.\n\n- The emissions are only from a single date. Emissions vary from day to day, and this information is lost. Similarly activity data is also not of a single day instead it is just a summation of all the activity from the power plants. \n"},{"metadata":{},"cell_type":"markdown","source":"### Feature Engineering\n\nTo assess what factors contribute to emissions, a model can be built to predict the NO2 emissions. The features of the model can be activities that sources of NO2, and/or causes of the spread of NO2. \n\nWe will create features from the NO2 emissions, weather, land assimilation data, and power plant data. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#Get the date information from each files name\ndef get_dates(file_path, data_source):\n    if data_source == 's5p':\n        fname_only = file_path.split('/')\n        dates_only = (fname_only[-1].split('_')[2], fname_only[-1].split('_')[3])\n        start_date = dates_only[0][:8]\n        end_date = dates_only[1][:8]\n        return start_date, end_date\n    elif data_source == 'gfs':\n        file_name = file_path.split('/')[-1]\n        date_w_extension = file_name.split('_')[-1]\n        date = date_w_extension.split('.')[0][:8]\n        return date\n    elif data_source == 'gldas':\n        file_name = file_path.split('/')[-1]\n        date = file_name.split('_')[1]\n        return date\n    \n#Read the data from the files and assemble a dataframe \ndef read_data(data_path,data_source):\n    data = []\n    for file in os.listdir(data_path):\n        file_path = data_path+file\n        img = tiff.imread(file_path)\n        img_to_add = {}\n        if data_source == 's5p':\n            start_date, end_date = get_dates(file_path, 's5p')\n            img_to_add['start_date'] = start_date\n            img_to_add['end_date'] = end_date\n            img_to_add['no2_emissions_mean'] = np.nanmean(img[:,:,0:4])\n        elif data_source == 'gfs':\n            img_to_add['date'] = get_dates(file_path, 'gfs')\n            img_to_add['temp'] = img[:,:,0]\n            img_to_add['specific_humidity'] = img[:,:,1]\n            img_to_add['relative_humidity'] = img[:, :, 2]\n            img_to_add['u_component_wind'] = img[:, :, 3]\n            img_to_add['v_component_wind'] = img[:, :, 4]\n            img_to_add['total_precipation'] = img[:, :, 5]\n        elif data_source == 'gldas':\n            img_to_add['date'] = get_dates(file_path, 'gldas')\n            for band in range(1, 13):\n                img_to_add['band'+str(band)] = rio.open(file_path).read(band)    \n                \n        data.append(img_to_add)\n        data_df = pd.DataFrame(data)\n        \n        if data_source == 's5p':\n            data_df['start_date'] = pd.to_datetime(data_df['start_date'])\n            data_df['end_date'] = pd.to_datetime(data_df['end_date'])\n            data_df.sort_values('start_date', inplace = True)\n            data_df.reset_index(drop = True, inplace = True)\n        else: \n            data_df['date'] = pd.to_datetime(data_df['date'])\n            data_df.sort_values('date', inplace=True)\n       \n    return data_df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"##### Getting Labels for Regression Model from Sentinel Satellite"},{"metadata":{"trusted":true},"cell_type":"code","source":"no2_path = eie_data_path+'/s5p_no2/'\nprint(\"Number of sentinel satellite pictures: \", len(os.listdir(no2_path)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below, we'll get the the data from .tif files that contain NO2 measurements of the Sentinel Satellite. We'll create a pandas dataframe that will have the mean no2 emissions. "},{"metadata":{"trusted":true},"cell_type":"code","source":"emissions_mean = read_data(no2_path, 's5p')\nemissions_mean.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below is a plot of the emissions mean (Y axis), and the start date on the x axis. We see that emissions peak just before May, followed by additonal peaks during June, July, and August."},{"metadata":{"trusted":true},"cell_type":"code","source":"emissions_mean.plot.line(x='start_date', y='no2_emissions_mean')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"##### Features from Atmosphere Data (Global Forecast System 384-Hour Predicted Atmosphere Data) "},{"metadata":{"trusted":true},"cell_type":"code","source":"gfs_path = eie_data_path + '/gfs/'\nprint(\"We have \" , len(gfs_path), \" pictures of the global forecast system\")\nweather_df = read_data(gfs_path, 'gfs')\nweather_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_weather_feature_stats(df):\n    for c in [col for col in df.columns if col not in ['date']]:\n        df[c] = df[c].apply(np.mean)\n    return df\nweather_mean_stats = get_weather_feature_stats(weather_df.copy())\nweather_mean_stats.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the DataFrame above you can see that each date as a atmosphere feature associated with it, i.e. temperature, humidity, etc. "},{"metadata":{},"cell_type":"markdown","source":"##### Beginning to assemble the dataset for our model\n\nBelow, we'll start building our dataset for our model. The dataset will have m examples in total, where m is the number of dates that correspond to the NO2 emissions measurements. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#As there are multiple weather measurements per day, group them by the date and take the mean\nweather_features = weather_mean_stats.groupby('date').mean()\n#Also group the emissions  by date \nemissions_mean = emissions_mean.groupby('start_date').mean()\n#concatenate the mean emissions dataframe and the weather features\ntraining_data = pd.concat([weather_features, emissions_mean], axis=1 , join='outer')\n#Drop any NaN values\ntraining_data.dropna(how='any',inplace=True)\ntraining_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"##### Features from Global Land Data Assimilation System"},{"metadata":{"trusted":true},"cell_type":"code","source":"gldas_path = eie_data_path + '/gldas/'\nprint(\"We have \", len(os.listdir(gldas_path)), \" pictures of land data\")\ngldas_df = read_data(gldas_path, 'gldas')\ngldas_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_gldas_mean_stats(df):\n    cols = [col for col in df.columns if col not in ['date']]\n    print(cols)\n    for c in cols:\n        df[c] = df[c].apply(np.mean)\n    return df\ngldas_mean_stats = get_gldas_mean_stats(gldas_df.copy())\ngldas_mean_stats['date'] = gldas_mean_stats['date'].apply(pd.Timestamp.date)\ngldas_mean_stats = gldas_mean_stats.fillna(0)\ngldas_mean_stats = gldas_mean_stats.groupby('date').mean()\n\ntraining_data = pd.concat([training_data, gldas_mean_stats], axis=1 , join='outer')\ntraining_data.dropna(how='any', inplace=True)\n#Add the total power generation of all the power plants\ntraining_data['total_power_generation'] = total_power_generation_2017\ntraining_data.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Linear Regression\n\nFor a regression model, our goal is to predict the NO2 emissions by using features engineered from power plant, atmosphere, and land use data. There are many factors that contribute to NO2 emissions. While there are sources of pollution such as electricty/power generation, traffic, and other industrial activities, conditions of land and the atmosphere can affect how pollutants remain in, and spread from the environment that they originate from.\n\nFor each given day we have the average of the NO2 emissions of that day, and the weather and land information as well. Features from weather data include the mean temperature, specific humidity, relative humidity, the U component of the wind, the V component of the wind, and the total precipation at surface. Similarly, the features from the land assimilation data include albedo, average surface skin temperature, plant canopy surface water, canpoy water evaporation, evaporation from bare soil, evapotranspiration, downward long-wave radiation flux, net long-wave radiation flux, potential evaporation rate, pressure, specific humidity, and heat flux (first 12 bands).There is not time series data available- so the total power generation of the entire island of Puerto Rico will be assumed to be constant (for the time being).\n\nMeasurements from s5p (NO2 emissions), weather, and land originally provide data for the entire island. For the purpose of the regression, we average the measurement accross the entire island. This does result in loss of spatial information, and is a much more coarse."},{"metadata":{"trusted":true},"cell_type":"code","source":"ds = training_data.copy()\ntrain_df = ds.iloc[:324]\ntest_df = ds[324:]\nX = train_df[[col for col in ds.columns if col not in ['no2_emissions_mean']]]\nY = train_df['no2_emissions_mean']\nmodel = LinearRegression().fit(X,Y)\ntest_X = test_df[[col for col in ds.columns if col not in ['no2_emissions_mean']]]\ntest_Y = test_df['no2_emissions_mean']\npreds = model.predict(test_X)\npreds = pd.Series(preds)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(preds.index, preds.values, color='red')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(test_Y.index, test_Y.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rms = np.sqrt(mean_squared_error(test_Y, preds))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(rms)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Based on the above we can see that there is room for improvement for the regression model. Future model improvements can include using more features from power plant/activity, weather, and land use data in addition to incorporating spatial information of emissions."},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}