{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction\nThis notebook introduces the statistical test performed on `omi` and how we constructed `df_ras_gen`."},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install rasterstats --quiet","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport glob \nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom datetime import datetime\n\nimport rasterio as rio\nimport folium\nimport tifffile\n\n# Geospatial libraries that we will be using for this\nimport rasterio\nimport rasterstats\nimport folium\nimport geopandas as gpd\nfrom shapely.geometry import Point\n\nfrom tqdm.notebook import trange, tqdm\n\nimport warnings\nwarnings.simplefilter(action='ignore')\n\npd.options.display.max_columns = None","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA_DIR = '../input/ds4g-environmental-insights-explorer/eie_data/'\nos.listdir(DATA_DIR)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We downloaded [EIA dataset](https://www.eia.gov/electricity/data/eia923/) for Puerto Rico."},{"metadata":{"trusted":true},"cell_type":"code","source":"eia_df_2018 = pd.read_excel('../input/eia-data/eia_puerto_rico_2018.xlsx')\neia_df_2017 = pd.read_excel('../input/eia-data/eia_puerto_rico_2017.xlsx')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We manually matched the plants that exist in EIA to `power_plants` from `gppd`. We only matched 22 out of 35 plants."},{"metadata":{"trusted":true},"cell_type":"code","source":"manual_map = {'AES ILUMINA': 'A.E.S. Corp.',\n                'AES Puerto Rico': 'AES Ilumina',\n                'Aguirre Plant': 'Aguirre',\n                'Cambalache Plant': 'Cambalache',\n                'Caonillas': 'Caonillas 1',\n                'Central San Juan Plant': 'San Juan CC', #'Caonillas 2'\n                'Costa Sur Plant': 'Costa Sur', # Carite 1\n                # 'Culebra': ''\n                'Daguao': 'Daguao',\n                'EcoElectrica': 'EcoEléctrica', # 'Dos Bocas'\n                # 'HIMA San Pablo - Caguas': ''\n                # 'Holsum de Puerto Rico, Inc.': 'Garzas 1'\n                # 'Humacao Solar Project, LLC': 'Garzas 2'\n                # 'Janssen Ortho LLC': 'Jobos'\n                'Jobos': 'Jobos', # Loiza Solar Park\n                'Mayaguez Plant': 'Mayagüez',\n                'Oriana Energy Hybrid': 'Oriana Solar Farm',\n                'Palo Seco Plant': 'Palo Seco',\n                'Pattern Santa Isabel LLC': 'Santa Isabel Wind Farm',\n                'Punta Lima Wind Farm': 'Punta Lima', # 'Río Blanco'\n                'San Fermin Solar Farm': 'San Fermin Solar Farm', # 'Salinas'\n                'Vega Baja': 'Vega Baja',\n                'Vieques': 'Vieques EPP', # 'San Juan CC'\n                'Yabucoa': 'Yabucoa'}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\npower_plants_path = '../input/gppd-geojson/Geolocated_gppd_120_pr.geojson'\npower_plants_original = gpd.read_file(power_plants_path)\n\n# we removed plants that do not exist in EIA data\npower_plants_path = '../input/geojson-matched/Geolocated_gppd_120_pr_matched.geojson'\npower_plants = gpd.read_file(power_plants_path)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# eia_data = pd.read_excel('../input/eiadata/eia_puerto_rico_2018.xlsx')\n\neia_df_2018 = pd.read_excel('../input/eia-data/eia_puerto_rico_2018.xlsx')\neia_df_2017 = pd.read_excel('../input/eia-data/eia_puerto_rico_2017.xlsx')\n\ndef process_eia_df(eia_data, power_plants):\n    eia_data.columns = eia_data.columns.str.replace('\\n', ' ')\n    \n    eia_data['mapped_name'] = eia_data['Plant Name'].map(manual_map)\n\n    eia_df = pd.merge(eia_data, power_plants, left_on='mapped_name', right_on='name', how='inner')\n\n    return eia_df\n\ndef calc_netgen(eia_df):\n    netgen_months = [col for col in eia_df.columns if col.find('Netgen') != -1]\n\n    netgen_months.extend(['system:index', 'primary_fuel'])\n\n    tmp = eia_df[netgen_months]\n\n    tmp = tmp.replace('.', 0)\n\n    netgen_sum_month = tmp.sum(axis=0)[:-2].values\n    \n    return netgen_sum_month\n    \ndef plot_netgen(eia_df_2017, eia_df_2018):\n    netgen_sum_month = []\n    netgen_sum_month.extend(calc_netgen(eia_df_2017))\n    netgen_sum_month.extend(calc_netgen(eia_df_2018))\n    \n    time = [str(i)+'-17' for i in range(1, 13)]\n    time.extend([str(i)+'-18' for i in range(1, 13)])\n    \n    # plt.figure(figsize=(8, 6))\n    plt.xticks(rotation=45)\n    plt.title('Total electricity generated from 2017 to 2018')\n    plt.xlabel('Time')\n    plt.ylabel('Electricity generated')\n    plt.plot(time, netgen_sum_month)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"eia_df_2018 = process_eia_df(eia_df_2018, power_plants)\neia_df_2017 = process_eia_df(eia_df_2017, power_plants)\n\nplot_netgen(eia_df_2017, eia_df_2018)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def filter_nan(eia_df):\n    df_sub = eia_df[eia_df['mapped_name'].isna()]\n    df_sub_filled = eia_df[~eia_df['mapped_name'].isna()]\n    return df_sub, df_sub_filled\n\ndf_sub, df_sub_filled = filter_nan(eia_df_2018)\n\ndf_sub_filled.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://disc.gsfc.nasa.gov/datasets/OMNO2_003/summary?keywords=OMI\nomi = pd.read_csv(\"../input/omiv1/OMI-Aura_L2-OMNO2.csv\")\n\nomi.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 1.0 Helper Functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"def scale_lat_to_img_idx(lat):\n    lat = float(lat)\n    lat_img_idx = (18.6 - lat) * 148 / (18.6 - 17.9)\n    return int(lat_img_idx)\n\ndef scale_lon_to_img_idx(lon):\n    lon = float(lon)\n    lon_img_idx = (67.3 + lon) * 475 / (67.3 - 65.2)\n    return int(lon_img_idx)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2.0 Validate External Data\nCheck if NO2 levels of +/- `DEFINED_BOUNDS` pixels from oil, gas, or coal plants are higher than other areas."},{"metadata":{"trusted":true},"cell_type":"code","source":"omi['img_idx_lt'] = omi['Latitude'].apply(scale_lat_to_img_idx)\n\nomi['img_idx_lg'] = omi['Longitude'].apply(scale_lon_to_img_idx)\n\nprint('len of omi before filter:', len(omi))\n\n# filter out of bounds\nomi = omi[(omi['img_idx_lt'] >= 0) & (omi['img_idx_lt'] <= 148)]\nomi = omi[(omi['img_idx_lg'] >= 0) & (omi['img_idx_lt'] <= 475)]\n\nomi.reset_index(inplace=True, drop=True)\n\nprint('len of omi after filter:', len(omi))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants['img_idx_lt'] = power_plants['latitude'].apply(scale_lat_to_img_idx)\n\npower_plants['img_idx_lg'] = power_plants['longitude'].apply(scale_lon_to_img_idx)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"non_rewew_energy = ['Oil', 'Coal', 'Gas']\n\nnon_renew_energy_plants = power_plants[power_plants['primary_fuel'].isin(non_rewew_energy)]\nrenew_energy_plants = power_plants[~power_plants['primary_fuel'].isin(non_rewew_energy)]\n\nnon_renew_energy_plants.head()\n\nlen(non_renew_energy_plants), len(renew_energy_plants)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How close are renewable energy plants to non-renewable energy plants?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot on the map\nlat = 18.200178; lon = -66.664513\n\ndef plot_polygons_on_Map(geodataframe,\n                         original_latitude=lat,\n                         original_longitude=lon,\n                         zoom=9):\n    plot = folium.Map(location = (original_latitude, original_longitude), zoom_start=zoom)\n\n    geojson = folium.GeoJson(geodataframe,\n                             style_function=lambda x: {\n                                 'Color':'white',\n                                 'fillColor':'red' if x['properties']['primary_fuel'] in ['Coal', 'Oil', 'Gas']\n                                                   else 'white'})\n    popup = folium.Popup(str(geodataframe.primary_fuel))\n    popup.add_to(geojson)\n    geojson.add_to(plot)\n\n    return(plot)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"buffered_power_plants = power_plants.copy()\n\nbuffered_power_plants['geometry'] = power_plants.geometry.buffer(0.1)\n\nplot_polygons_on_Map(buffered_power_plants)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We mark an area as a `is_polluted_bounds` area when the pixel is +/- 15 pixels from a non renewable energy plant."},{"metadata":{"trusted":true},"cell_type":"code","source":"polluted_idx = non_renew_energy_plants[['img_idx_lt', 'img_idx_lg']].values\n\nnon_polluted_idx = renew_energy_plants[['img_idx_lt', 'img_idx_lg']].values\n\nexamples = omi[['img_idx_lt', 'img_idx_lg']].values\n\nDEFINED_BOUNDS = 15\n\ndef in_polluted_bounds(polluted_idx, lat, lon):\n    for pol_lat, pol_lon in polluted_idx:\n        lat_in_bounds = pol_lat - DEFINED_BOUNDS <= lat <= pol_lat + DEFINED_BOUNDS\n        lon_in_bounds = pol_lon - DEFINED_BOUNDS <= lon <= pol_lon + DEFINED_BOUNDS\n        if lat_in_bounds and lon_in_bounds:\n            return True\n    return False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# omi['is_polluted_bounds'] = omi.apply(lambda x: in_polluted_bounds(polluted_idx, x['img_idx_lat'], x['img_idx_lon']), axis=1)\n\npolluted_bounds_bool = [0] * len(omi)\n\nfor idx, row_data in omi.iterrows():\n    lat, lon = row_data['img_idx_lt'], row_data['img_idx_lg']\n    if in_polluted_bounds(polluted_idx, lat, lon):\n        polluted_bounds_bool[idx] = 1\n\nomi['is_polluted_bounds'] = polluted_bounds_bool\n\nnon_polluted_bounds_bool = [0] * len(omi)\n\nfor idx, row_data in omi.iterrows():\n    lat, lon = row_data['img_idx_lt'], row_data['img_idx_lg']\n    if in_polluted_bounds(non_polluted_idx, lat, lon):\n        non_polluted_bounds_bool[idx] = 1\n\nomi['not_polluted_bounds'] = non_polluted_bounds_bool","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# omi.dropna(inplace=True)\n\nomi['is_polluted_bounds'].dropna(inplace=True)\nomi['not_polluted_bounds'].dropna(inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 1. Statistical Test for `omi`"},{"metadata":{"trusted":true},"cell_type":"code","source":"omi.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import stats\n\ndef run_t_test(population_1, population_2, num_tests=10, num_samples=500):\n    count = 0\n\n    for i in range(num_tests):\n        np.random.seed(i)        \n\n        sample_polluted = np.random.choice(population_1, num_samples)\n        sample_non_polluted = np.random.choice(population_2, num_samples)\n\n        _, p_value = stats.ttest_ind(sample_polluted, sample_non_polluted)\n\n        if p_value < 0.05:\n            count += 1\n            \n    return count / num_tests","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"col = 'ColumnAmountNO2'\n\ndirty_no2 = omi[omi['is_polluted_bounds'] == 1][col].mean() \nclean_no2 = omi[omi['is_polluted_bounds'] == 0][col].mean()\n\npolluted_higher = (dirty_no2 - clean_no2) / clean_no2\npct_stats_signific = run_t_test(omi[omi['is_polluted_bounds'] == 1][col], omi[omi['is_polluted_bounds'] == 0][col])\n\npolluted_higher, pct_stats_signific","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dirty_no2 = omi[omi['is_polluted_bounds'] == 1][col].mean() \nclean_no2 = omi[omi['not_polluted_bounds'] == 1][col].mean()\n\npolluted_higher = (dirty_no2 - clean_no2) / clean_no2\npct_stats_signific = run_t_test(omi[omi['is_polluted_bounds'] == 1][col], omi[omi['not_polluted_bounds'] == 1][col])\n\npolluted_higher, pct_stats_signific","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* NO2 at area +/- 15 pixels from a `Coal`, `Gas` or `Oil` plant is 6.96% higher than other regions. \n* NO2 at area +/- 15 pixels from a `Coal`, `Gas` or `Oil` plant is 0.6% higher than a plant that is +/- 15 pixels from a `Hydro`, `Solar`, and `Wind` plant.\n\nWith a sample of 500 data points, the first result is statistically significant for 8 out of 10 random sampling while the second result is **not** statistically significant."},{"metadata":{},"cell_type":"markdown","source":"## 2. Constructing `df_ras_gen`"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_s5 = pd.DataFrame()\n\nfiles = []\ncaptured_datetime = []\n\nfor dirname, _, filenames in os.walk(os.path.join(DATA_DIR, 's5p_no2')):\n    for filename in filenames:\n        files.append(os.path.join(dirname, filename))\n        captured_datetime.append(datetime.strptime(filename.split('_')[-2],'%Y%m%dT%H%M%S').date())\n        \ndf_s5['path'] = files\ndf_s5['captured_date'] = captured_datetime\n\n# Sort dataframe by ascending time\ndf_s5.sort_values('captured_date', inplace=True)\ndf_s5.reset_index(inplace=True, drop=True)\n\ndf_s5.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# https://www.kaggle.com/maxlenormand/saving-the-power-plants-csv-to-geojson\npower_plants_path","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_df(df_s5, power_plants_path, df, BAND=1):\n    for idx, row_data in tqdm(df_s5.iterrows()):\n        raster_stats = rasterstats.zonal_stats(power_plants_path,\n                                               row_data['path'],\n                                               band=BAND,\n                                               stats=['mean'])\n\n        raster_stats_list = [plant['mean'] for plant in raster_stats]\n\n        df[row_data['captured_date']] = raster_stats_list\n        \n    return df\n\n# N02_measurments_df = pd.DataFrame({'index_power_plant': power_plants['system:index']})\n# N02_measurments_df = create_df(df_s5, power_plants_path, N02_measurments_df, 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We save each of the band into a `csv` file and reload the files under `rasterstats`."},{"metadata":{"trusted":true},"cell_type":"code","source":"features = ['NO2_column_number_density',\n            'tropospheric_NO2_column_number_density',\n            'stratospheric_NO2_column_number_density',\n            'NO2_slant_column_number_density',\n            'tropopause_pressure',\n            'absorbing_aerosol_index',\n            'cloud_fraction',\n           ]\n\ntmp = pd.DataFrame({'index_power_plant': power_plants['system:index']})\n\nfor idx, feat in enumerate(features):\n    df = create_df(df_s5, power_plants_path, tmp, idx+1)\n    \n    filename = feat + '.csv'\n    \n    df.to_csv(filename, index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def calc_avg_month(df, feature):\n    power_plants = df.iloc[:, 0]\n    tmp = df.iloc[:, 1:]\n    tmp.columns = pd.to_datetime(tmp.columns)\n    \n    # transpose & reset\n    tmp = tmp.T\n    tmp.reset_index(inplace=True)\n    \n    # extract year\n    tmp['year'] = tmp['index'].dt.year\n    tmp['month'] = tmp['index'].dt.month\n    \n    df = pd.DataFrame()\n    \n    dict_name = {\n        'min': feature + '_min',\n        'max': feature + '_max',\n        'mean': feature + '_mean',\n        'std': feature + '_std'\n    }        \n            \n    for idx, plant in enumerate(power_plants):\n        subset = tmp.groupby(['year', 'month'])[idx].agg(['min', 'max', 'mean', 'std'])\n        subset = subset.reset_index()\n        subset['system:index'] = plant   \n        subset = subset.rename(columns=dict_name)\n        df = pd.concat([df, subset], axis=0, ignore_index=True)\n            \n    return df","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here, we reload the files for each band from Sentinel-5P OFFL NO2."},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_and_merge():\n    raster_files_dir = '../input/rasterstats/rasterstats'\n\n    raster_files = os.listdir(raster_files_dir)\n    \n    comb_df = pd.DataFrame()\n    \n    for file in raster_files:\n        df = pd.read_csv(os.path.join(raster_files_dir, file))\n        \n        feature, _ = file.split('.')\n        \n        df = calc_avg_month(df, feature)\n        \n        comb_df = pd.concat([comb_df, df], axis=1)\n            \n    return comb_df.T.drop_duplicates().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_ras = read_and_merge()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Reorder columns."},{"metadata":{"trusted":true},"cell_type":"code","source":"cols = df_ras.columns\n\nnew_cols = ['system:index']\n\nnew_cols.extend([c for c in cols if c != 'system:index'])\n\ndf_ras = df_ras[new_cols]\n\ndf_ras.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Merge EIA data with `df_ras`\n\nFrom the EIA data, we have the monthly electricity generated for each plant. We then merge it with `df_ras` that contains the `mean`, `max`, `min`, and `std` rasterstats."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_ras['generated'] = np.nan\n\ndf_ras.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def safe_convert(x):\n    try:\n        x = float(x)\n        return x\n    except:\n        return np.nan\n\n    \ndf_sub, df_sub_filled = filter_nan(eia_df_2018)\n\nnetgen_months = [col for col in df_sub_filled.columns if col.find('Netgen') != -1]\n\ntmp = df_ras.copy()\n\nfinal = pd.DataFrame()\n\nfor idx, m in enumerate(netgen_months[6:], 6):\n    df_sub_filled[m] = df_sub_filled[m].apply(safe_convert)\n    \n    g = df_sub_filled.groupby(['system:index'])[m].sum()\n    \n    g = g.reset_index()\n    \n    for _, row in g.iterrows():\n        name_con = tmp['system:index'] == row['system:index']\n        month_con = tmp['month'] == idx + 1\n        tmp.loc[(name_con) & (month_con), 'generated'] = row[m]\n    \ntmp.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_sub, df_sub_filled = filter_nan(eia_df_2017)\n\nnetgen_months = [col for col in df_sub_filled.columns if col.find('Netgen') != -1]\n\nfor idx, m in enumerate(netgen_months[:6]):\n    df_sub_filled[m] = df_sub_filled[m].apply(safe_convert)\n    \n    g = df_sub_filled.groupby(['system:index'])[m].sum()\n    \n    g = g.reset_index()\n    \n    for _, row in g.iterrows():\n        name_con = tmp['system:index'] == row['system:index']\n        month_con = tmp['month'] == idx + 1\n        tmp.loc[(name_con) & (month_con), 'generated'] = row[m]\n        \ntmp.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tmp.to_csv('df_ras_gen.csv', index=False)","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}