{"cells":[{"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\nfrom sklearn.preprocessing import PolynomialFeatures","execution_count":null,"outputs":[]},{"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":{"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":{"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":{"trusted":true},"cell_type":"code","source":"print('There are {} power plants'.format(gpp_df.shape[0]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"จำนวนโรงงานทั้งหมดมี 35 โรงงาน"},{"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":"เชื้อเพลิงไหนที่โรงงานสร้างเยอะที่สุด"},{"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":"หาว่าเชื้อเพลินไหนผลิตไฟฟ้าได้เยอะที่สุดในเกาะ"},{"metadata":{"trusted":true},"cell_type":"code","source":"gpp_df['commissioning_year'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"เช็คอายุโรงงาน"},{"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":"หาค่าควจุของเชื้อเพลิงแต่ชนิด"},{"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":{"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":{"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":{"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":{"trusted":true},"cell_type":"code","source":"tiff.imread(image).shape","execution_count":null,"outputs":[]},{"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":"ในgoogle ได้บอกไว้ว่า ค่าNO2 ที่ปล่อยจากโรงงานนั้นคิดเป็นแค่ 14% ของNO2 ทั้งหมด\nE / A = EF\n\nE = emissions\nA = activity rate\nEF = emission factor"},{"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":"มีภาพทั้งหมด 387 ภาพ\n"},{"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":{"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":{"trusted":true},"cell_type":"code","source":"Check จำนวน Nan Value ทั้งหมดที่มีในแต่ละ Band และทำการ ignored","execution_count":null,"outputs":[]},{"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":{"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":"หลังจาก ignored แล้วเอา มา plot ค้นพบว่าบางเวลามีค่า NO2 ต่ำเกินเพราะว่าเรามีข้อมูลไม่ครบของวันนั้นๆ"},{"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":"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":{"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":"มีข้อมูลซ้ำกัน 42 วัน"},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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":"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 / p_type_sum\nprint(f'Simplified emissions factor for Puerto Rico is {factor} mol * h / m^2 * gw')","execution_count":null,"outputs":[]},{"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":{"trusted":true},"cell_type":"code","source":"tiff.imread(weather_pictures_path[0]).shape","execution_count":null,"outputs":[]},{"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":{"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":{"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":{"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":"หาค่า colleration ออกมาเป็น Heatmap สังเกตุว่าค่าที่สัมพันธ์กับ NO2 มากที่สุดคือ Temperature กับ Humidity เลยเลือก 2 Attribute นี้ไปทำการ Train model"},{"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', \n                                '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\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    error = 0\n    poly_features = PolynomialFeatures(2, include_bias = True)\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        print(np.shape(x_train))\n        print(np.array(x_train))\n        print(np.shape(x_val))\n        print(np.array(x_val))\n        print(np.shape(y_train))\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        print(np.shape(x_train))\n        print(np.shape(x_val))\n        #Polynomial Feature\n        print(np.shape(x_train))\n        print(np.shape(x_val))\n        x_poly = poly_features.fit_transform(x_train)\n        print(np.shape(x_poly))\n        print(np.array(x_poly))\n        x_val = poly_features.transform(x_val)\n        print(np.shape(x_val))\n        test_scaled = poly_features.transform(test[features])\n        month = train[['month']].iloc[val_ind]['month'].unique()[0]\n        model = LinearRegression().fit(x_poly, y_train)\n        r2=model.score(x_poly,y_train)\n        fold_prediction = model.predict(x_val)+(error*0.1)\n        print(np.shape(fold_prediction))\n        print(np.shape(oof[val_ind]))\n        fold_error = np.sqrt(metrics.mean_squared_error(y_val, fold_prediction))\n        error = y_val - fold_prediction\n        error = np.mean(error)\n        print(\"ssssss\" ,r2)\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)\n#plot_coef(all_coef1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"ในการทำ Train model ในนี้ใช้ Polynomial Regression โดย Set degree เป็น 2 และ K fold = 12 ครั้ง แล้วทำการหา RMSE ของ Train และ Test ซึ่งได้เท่ากับ\nOur out of folds rmse is 1.55\nOur test rmse is 1.79"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""}],"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}