{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Methodology for the Emissions Factors Calculation combining a physical approach with satellite images\n\n\nIn this notebook, we suggest a physical approach to compute the emission factor per power plant and demonstrate this on the example of Puerto Rico dataset. This approach is based on a physical model to simulate the spreading of NO2 in the atmosphere due to the diffusion and the convection. By using the superposition principle, we can decompose the solution to a linear combination of a set of solutions that do not depend on the emission rate of each power plant. By using this trick, we can then find the best emission rate for each power plant by solving a simple constrained optimization problem. The optimization consists in finding the best fit of emission rate parameters with respect to satellite images. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Plan:\n1. <u>**Data analysis**</u>\n   1. **Puerto rico power plant:**\n       \n   2.  **S5P No2 analysis**\n       * correlation between No2 concentration and the monthly generated electricity \n       * No2 concentration and power plant location\n   3. **GFS data (wind) analysis**\n       * Wind impact on No2 concentration\n \n       \n2. <u>**Modeling:**</u>\n\n    1. ** Modeling of No2 spreading in the atmosphere**\n        - Diffusion-convection with source term of No2\n        - Turbulence modeling\n    2. **Modeling of the emission rate:**\n       - Basis decomposition of No2 concentration\n       - Least square fit modeling of the emission rate\n   \n3. <u>**Validation:**</u>\n  \n  to be done","execution_count":null},{"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\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport sys\nimport os\nfrom matplotlib import animation, rc\nimport datetime\nimport rasterio as rio\nimport folium\n#from scipy.ndimage import gaussian_filter\nfrom utils_ds4g import *\nfrom no2_solver import *\nsys.path.append('')\nmain_folder = '/kaggle/input/ds4g-environmental-insights-explorer/'\nfolder_s5p_no2 = main_folder+'eie_data/s5p_no2/'\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 os\nimport sys\nsys.path.append('/kaggle/input/ds4g-environmental-insights-explorer/')\n\n\n\n\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Important parameter\nshape0 = (148, 475)\nlat_limits = np.array([ 17.9, 18.56])\nlon_limits = np.array([-67.32, -65.19])\n[lat, lon] = [lat_limits.mean(),  lon_limits.mean()]\nRt = 6371*1.e3\nT_year_h = 24*365.25\nT_day_s = 24*3600\n\nT_year_s = T_year_h*3600\n# convert the polor coordiate to a local cartesian coordinate by assuming a locally flat area \nX = Rt*np.cos(lat*np.pi/180)*np.linspace(np.pi/180*lon_limits[0],np.pi/180*lon_limits[1],shape0[1])\nY = Rt*np.linspace(np.pi/180*lat_limits[0],np.pi/180*lat_limits[1],shape0[0])\ndx = np.abs(X[1]-X[0])\ndy = np.abs(Y[1]-Y[0])\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 1- Data analysis\n## 1.1 Existing power plant in Puerto rico","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"power_plants = pd.read_csv(main_folder+'eie_data/gppd/gppd_120_pr.csv')\npower_plants.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* **Generated energy per year [2013-2017]** columns are not given and only the estimated one is reported\n* <b>fuel_1, fuel_2 and fuel_3</b> are not used in all powerplants and only primary_fuel is relevant here\n* The only considered country is US, therefore we do not need to keep **country data** \n\n-> Let's do some cleaning","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# drop useless columns\nfor y in range(2013,2018):\n    del power_plants['generation_gwh_'+str(y)]\nfor i in range(1,4):\n    del power_plants['other_fuel'+str(i)]\ndel power_plants['country']\ndel power_plants['country_long']\npower_plants.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# fuel mix \nprint(power_plants.primary_fuel.value_counts())\nsns.countplot(power_plants.primary_fuel)\n\n# the total capacity of power plant by fuel type\n\n\npower_plants.groupby('primary_fuel',group_keys=True).sum().reset_index().plot.bar(x='primary_fuel',y='capacity_mw')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Puerto rico has a mixed energy system. However, thermal power plants, which use fossil fuel (oil, coal and gas), have dominant capacity compared to \"green power plant\" (Hydro, wind and solar). Since our interest in this study is the emission factor, we will consider in the following only thermal power plants.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"thermal_power_plants = power_plants[power_plants.primary_fuel.isin(['Coal','Gas','Oil'])].reset_index()\nprint('percentage of fossil fuel power plant capacity',thermal_power_plants.capacity_mw.sum()/power_plants.capacity_mw.sum())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# capacity factor = total generated energy in one year/ (maximal energy that can be generated over the year)\n\nthermal_power_plants['capacity_factor'] = 1.e3*thermal_power_plants['estimated_generation_gwh']/ (T_year_h*thermal_power_plants['capacity_mw'])\nthermal_power_plants.plot.bar(x='primary_fuel',y='capacity_factor')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### The capacity factor is defined as the ratio of generated electricity in a given period divided by the maximal enery that we can generate in the same periode. Therefore, capacity facot is included between 0 and 1.0. The capacity factor calculated for the coal power plant gives a very high value more than 100. We need to correct the emission factor by correcting the estimated_generation_gwh (the capacity value of the coal looks Ok). The annual capacity factor of Coal power plant in USA (based on data in https://www.eia.gov/electricity/monthly/epm_table_grapher.php?t=epmt_6_07_a) is about 50%. In the following, we use this value to correct the estimated_generation_gwh.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#thermal_power_plants.loc[thermal_power_plants.primary_fuel=='Coal'].capacity_factor  = 0.5\nthermal_power_plants.at[14,'capacity_factor'] = 0.5 # we update the capacity factor of the coal power plant\nthermal_power_plants['estimated_generation_gwh'] =  (T_year_h*thermal_power_plants['capacity_mw'])*thermal_power_plants['capacity_factor']/(1.e3) \nthermal_power_plants.plot.bar(x='primary_fuel',y='capacity_factor')\n\nthermal_power_plants.groupby('primary_fuel',group_keys=True).sum().reset_index().plot.bar(x='primary_fuel',y='estimated_generation_gwh',title='total generated electricty per fuel')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Now, the capacity factor and total generated electricity looks fine.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"thermal_power_plants = split_column_into_new_columns(thermal_power_plants,'.geo','latitude',50,66)\nthermal_power_plants = split_column_into_new_columns(thermal_power_plants,'.geo','longitude',31,48)\nthermal_power_plants['latitude'] = thermal_power_plants['latitude'].astype(float)\na = np.array(thermal_power_plants['latitude'].values.tolist()) # 18 instead of 8\nthermal_power_plants['latitude'] = np.where(a < 10, a+10, a).tolist() \n\n\nplot_points_on_map(thermal_power_plants,0,425,'latitude',lat,'longitude',lon,9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Puerto ricco power plant analysis\n1. Puerto ricco has a mixed energy system, where the fossil fuel power plants consist about 94% of the total capacity, and this percentage will be higher for generated electricty (because \"green\" power plant suffer from low capacity factor due to intermittence).\n2. The fossil power plant are located along the coast where we expect strong wind flow. The wind will be one of the key parameters that need to be explored and taken in the model to estimate the average historical emission factor and the emission factor per power plant,\n3. From the data base of the thermal power plant, we will use in our study the capacity_mw, yearly estimated_energy_mw and location of the power plant. Another data that seems to me important is also the montly estimated_energy. Indeed, we want to study the correlation between the wind and NO2 concentration, this study can be misleading if we do not take into account the generated electricty in a smaller period than the whole year. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"###  Additional data: monthly generated electricity\n\nIn the following, we are going to add an estimation of the generated electricity each month. To estimate the generated electricity, we use data of monthly generated electrity in USA between 2018 and 2019(<url> https://www.eia.gov/electricity/data/browser/#/topic/0?agg=2,0,1&fuel=vtvv&geo=g&sec=g&freq=M&start=201801&end=201912&ctype=linechart&ltype=pin&rtype=s&maptype=0&rse=0&pin=</url>). We the compute the percentage of generated electrity per month with respect to the whole year. I will use this data to compute an estimation of the generated electicity per month. \n\nIn a second part of this section, an estimation of a monthly NO2 emission will be given based on the monthly estimation of the generated electrity and a generic emission factor of fossil power plant (","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# the monthly capacity factor are in https://www.eia.gov/electricity/data/browser/#/topic/0?agg=2,0,1&fuel=vtvv&geo=g&sec=g&freq=M&start=201801&end=201912&ctype=linechart&ltype=pin&rtype=s&maptype=0&rse=0&pin=\nmonthly_generated_electricity_USA = pd.read_csv('/kaggle/input/monthly-generated-electricity/Net_generation_for_all_sectors.csv')\nmonthly_generated_electricity_USA.set_index('description',inplace=True)\nmonthly_generated_electricity_USA = monthly_generated_electricity_USA.transpose().reset_index()\nmonthly_generated_electricity_USA = monthly_generated_electricity_USA.rename(columns={\"index\": \"date\"})\nmonthly_generated_electricity_USA.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"month_num_map = {'Jan':1,'Feb':2,'Mar':3,'Apr':4,'May':5,'Jun':6,'Jul':7,'Aug':8,'Sep':9,'Oct':10,'Nov':11,'Dec':12 }\n\nmonthly_generated_electricity_USA['month'] = monthly_generated_electricity_USA['date'].apply(lambda r: r.split(' ')[0])\nmonthly_generated_electricity_USA['month_num'] = monthly_generated_electricity_USA['month'].apply(lambda r: month_num_map[r])\n\nu = monthly_generated_electricity_USA.groupby('month',group_keys=True).mean()\nu.sort_values(by='month_num',inplace=True)\nu = u.reset_index()\n\nu['all_fuels_month_percent'] = u['all fuels']/(u['all fuels'].sum())\nfig,ax = plt.subplots(1,2,figsize=(14,6))\nu.drop(['month_num'],axis=1)\nu.plot(x='month',y='all fuels',ax=ax[0])\nax[0].set_ylabel('monthly generated electricity (gwh)')\nu.plot(x='month',y='all_fuels_month_percent',ax=ax[1])\nax[1].set_ylabel('month fraction of generated electricity')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### It will be for sure better if we can use monthly data of Puerto rico rather than of all USA. But, we can expected that the trend of the monthly variation of the generated electricity will be close to the one of Puerto-ricco. We will also suppose that the month fraction of generated electricity per power plant is as the global one.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\nfor [month, fraction] in zip(u.month,u.all_fuels_month_percent):\n    thermal_power_plants['estimated_generation_gwh_'+month] = thermal_power_plants['estimated_generation_gwh']*fraction\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The average emission factor of NO2 per fossil power plant can be found in this report ():\n* Coal: 1.6 g/MBTU\n* Oil: 0.6  g/MBTU\n* Gas: 0.1 g/MBTU\n\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"molar_mass_NO2 = 46.005 # g/mol\nMBTU_to_KWH = 0.29307\nEmission_factor_g_per_MBTU = {'Coal':1.6,'Gas':0.1,'Oil':0.6}\nEmission_factor_mol_per_kwh = Emission_factor_g_per_MBTU\nfor k in Emission_factor_g_per_MBTU.keys():\n    Emission_factor_mol_per_kwh[k] = Emission_factor_g_per_MBTU[k]/molar_mass_NO2*MBTU_to_KWH\n    \nthermal_power_plants['estimated_emission_factor_mol_kwh'] = thermal_power_plants.primary_fuel.apply(lambda r:Emission_factor_mol_per_kwh[r])\nfor [m, p] in zip(u.month,u.all_fuels_month_percent):\n    thermal_power_plants['estimated_emission_mol_'+m] = 1.e6*thermal_power_plants['estimated_emission_factor_mol_kwh']*thermal_power_plants['estimated_generation_gwh_'+m]\n\n    \n#  Plot the no2 emitted quantity per power plant in mol\nfig = plt.figure(figsize=(10,10))\ncolor_perfuel = {'Oil':'b','Coal':'k','Gas':'g'}\nfor i in range(len(thermal_power_plants)):\n    #fuel.append(thermal_power_plants.loc[i,'estimated_generation_gwh_Jan'])\n    emission_per_month = np.zeros(12)\n    for im in range(12):\n        emission_per_month[im] = thermal_power_plants.loc[i,'estimated_emission_mol_'+months_list[im]]\n    fuel_name = thermal_power_plants.loc[i,'primary_fuel']   \n    plt.plot(months_list,emission_per_month,color=color_perfuel[fuel_name],label=fuel_name)\nplt.legend()\nplt.xlabel('month')   \nplt.ylabel('No2 emission (mol)')  ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The current estimation of the monthly emission per power plant is a very rough estimation and we do not aim to use it as ground truth in the following, but rather as result to compare with the model that we are going to design in the following. This current estimation gives us an idea about the expected range of the emission values. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 1.2 - S5P No2 concentration","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"image = main_folder+'eie_data/s5p_no2/s5p_no2_20180708T172237_20180714T190743.tif'\nimage_band = rio.open(image)\nimage_band.descriptions","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sentinel 5P Offline NO2 has 12 channels (bands). The four first bands are the concentration in mol/m^2 in four different vertical columns. \n\nIn the following, we use the band 2  which corresponds to tropospheric_NO2_column_number_density. The tropospheric layer corresponds to the low layer of the atmosphere (the average height of the troposphere is 18 km), which is more relevant for the current study.\n\n\nRemark: the tiff images of s5p are flipped in Y-axis, you can refer to this notebook (https://www.kaggle.com/essadkim/are-tiff-images-flipped-in-y-axis).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"s5p_filenames =  np.sort(os.listdir(folder_s5p_no2))\nimages_s5p = {b:np.array([rio.open(folder_s5p_no2+file).read(b)[::-1,:] for file in s5p_filenames]) for b in range(1,5) }\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"start_date = datetime.date(2018, 7, 7)\nend_date = start_date+  datetime.timedelta(days=5)\n# we apply a gaussian filter and average over the periode start_date -> end_date\nim_no2 = nan_filter_mean(folder_s5p_no2,images_s5p,start_date,end_date,band=2,sigma=5.,acceptable_nan=0.1)\nfig = plt.figure(figsize=(8,8))\nplt.imshow(im_no2,cmap='jet')\nplt.colorbar()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_concentration_month = np.zeros(12)\ngenerated_energy = np.zeros(12)\ndate = []\ndebut_date_list = [  datetime.date(2018, i, 1) if (i<13) else datetime.date(2019, i-12, 1)   for i in range(7,20)]\n\nfor i in range(0,12):\n  \n    debut_date = debut_date_list[i]\n    end_date = debut_date_list[i+1]\n   \n    \n    im_no2 = nan_filter_mean(folder_s5p_no2,images_s5p,debut_date,end_date,band=2,sigma=0.,acceptable_nan=0.1)\n    month = debut_date.ctime().split(' ')[1]\n    date.append(month)\n    generated_energy[i]  = thermal_power_plants['estimated_generation_gwh_'+month].sum()\n    NO2_concentration_month[i] = im_no2.sum()\n    \nNO2_concentration_month_normalize = NO2_concentration_month/NO2_concentration_month.mean()\ngenerated_energy_normalize = generated_energy/generated_energy.mean()\n\nplt.plot(date,NO2_concentration_month_normalize,date,generated_energy_normalize)\nplt.legend(['No2 concentration','generated electricity'])\nprint('the correlation between NO2 concentration and the generated electricity is ',np.corrcoef(NO2_concentration_month,generated_energy)[0,1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- We can see that the concentration of NO2 is moderately correlated to the generated electricity with a factor  of 0.43\n- Other factors can have also an important impact on NO2 concentration:\n     + Other emission source (cars, aeroport,..)  \n     + Wind intensity and direction\n     + Temperature and humidty : NO2 can have some reactions with other chemical components in air like (H2o,O3,...) and their equilibrium reaction factor will depend on these physical factors \n        ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Power plant location and NO2 concentration","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\n# the mean no2 concentration on one year\n# the mean no2 concentration on one year\nstart_date = datetime.date(2018, 7, 7) \nend_date = start_date+  datetime.timedelta(days=365)\nmean_im_no2_year = nan_filter_mean(folder_s5p_no2,images_s5p,debut_date,end_date,band=2,sigma=5.,acceptable_nan=0.1)\n# load the power plant location and the generated elecricity on one year\n\nxy_location, elect_year,emission_per_month = power_plant_location(thermal_power_plants,mean_im_no2_year.shape,lat_limits,lon_limits)\n\nelect_year_norm = elect_year/elect_year.mean()\n\ns = mean_im_no2_year.shape\nfig = plt.figure( figsize=(12,9) )\nplt.imshow(mean_im_no2_year,cmap='jet')  # +0.2*np.mean(im_no2)*im_plants\nfor [xy,c] in zip(xy_location,elect_year_norm):\n    plt.plot(xy[0],xy[1],'o',markersize=5*c,color='k')\n    #plt.annotate('o', xy=xy, xycoords='data'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see a high concentration of No2 located in the center-north of Puerto ricco, where mainly two Oil power plants with high capacity are located: \n  -  Palo Seco,\tOil power plant of capacity 602MW,\n  -  San Juan CC,\tOil power plant of capacity 864.MW.\t \n  \nIn center-south, we can also see a moderate concetration of No2, where we have also two big power plant in terms of capacity:\n  - Aguirre, an Oil power plant of capacity 1492MW\n  - A.E.S. Corp., a Coal power plant of capacity 454MW.\n\nThe concentration in south is still very low compared to the north, even if the total capacity of the power plants in the two region are almost the same. More than that, in south, one of the power plant is a coal power plant, which its emission factor of NO2 is about three times more than Oil power plant. This emphasis that we need to take other physical parameter into account for reliable modeling of the emission No2 quantity per power plant.\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 1.3- GFS data\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"folder_gfs = main_folder+'eie_data/gfs/'\ngfs_filenames =  np.sort(os.listdir(folder_gfs))\nimages_gfs = {b:np.array([np.array(rio.open(folder_gfs+file).read(b)[::-1,:]) for file in gfs_filenames]) for b in range(4,6) }","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nb_days = 2. # we average data over nb_days of days\nnb_months = 10\nsigma = 5. # sigma to be applied for gaussian filter \nfig, ax = plt.subplots(nb_months,2, figsize=(60, 60))\n\nvmin, vmax = np.nanpercentile(mean_im_no2_year, (5,95))\n[vmin, vmax ] = [0.8*vmin,1.2*vmax]\n#[vmin, vmax ] = [-5.0*1.e-6,5.0e-6]\n\nfor i in range(0,nb_months):\n    debut_date = datetime.date(2018, 7, 7) + datetime.timedelta(days=i*30)\n    end_date = debut_date+  datetime.timedelta(days=nb_days)\n    no2f_filtered = nan_filter_mean(folder_s5p_no2,images_s5p,debut_date,end_date,band=2,sigma=sigma)\n    ax[i,0].imshow(no2f_filtered,cmap='jet',vmin=vmin,vmax=vmax)\n    \n   # ax[i,0].plot(xy[0],xy[1],'o',markersize=20*c,markerfacecolor=\"None\",\n    #     markeredgecolor='red', markeredgewidth=5)\n    \n   # im_prec = no2f_filtered\n    ax[i,0].set(title=str(debut_date)+' to '+str(end_date))\n    for [xy,c] in zip(xy_location,elect_year_norm):\n        ax[i,0].plot(xy[0],xy[1],'o',markersize=20*c,markerfacecolor=\"None\",\n         markeredgecolor='red', markeredgewidth=5)\n    #plt.colorbar()\n    # your code here    \n    \n    \n    \n    mean_u = nan_filter_mean(folder_gfs,images_gfs,debut_date,end_date,band=4,sigma=0.)\n    mean_v = nan_filter_mean(folder_gfs,images_gfs,debut_date,end_date,band=5,sigma=0.)\n    ## reduce the number of dimension \n    [h_step, w_step] = [20,40]\n    shape0 = mean_u.shape\n\n    mean_u_red = mean_u[::h_step,::w_step]\n    mean_v_red = mean_v[::h_step,::w_step]\n    shape = mean_u_red.shape\n    \n    x = np.linspace(lon_limits[0],lon_limits[1],shape0[1])\n    y = np.linspace(lat_limits[0],lat_limits[1],shape0[0])\n    ax[i,1].quiver(x[::w_step],y[::h_step],mean_u_red,mean_v_red,scale=40) \n    for [xy,c] in zip(xy_location,elect_year_norm):\n        col = x[int(xy[0])]\n        row = y[int(xy[1])]\n        ax[i,1].plot(col, row,'o',markersize=20*c,color='r')\n    ax[i,1].set(title=str(debut_date)+' to '+str(end_date))\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Wind analysis:\n\n\nThe main direction of the wind is from East to West. In some period of the year, we also see a part of the wind coming from the south-East than it is redirected (due to wind-resinstance of the land topographie) to the west. In the noth part of Puerto ricco, the wind stream is mainly from East to West over the whole year, with a small component of the wind speed toward the South (outside of the land). This can explain the low NO2 concentration in the north part pf the region. \n\nAnother important remark is that: NO2 concentration peaks are presence close to some power plants in the south region where the wind speed decelerate. The deceleration of the wind can cause an accumalation of aerosols drifted by the wind.\n\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"1. ### Correlation between wind, estimated generated electricity and no2 concentration\n\nThe correlation between No2 concentration in Puerto ricco and the monthly generated electricity is about 0.42 as was shown previousely. \nWe explain that this correlation can be high if we include the wind parameter. Indeed, in a windy day we expect a decrease of the No2 for the same produced electrity and an increase of No2 concentration in a non-windy day. \n\nTo explore the correlation bewteen generated electricity and No2 concentration, we propose to use the estimated monthly emission, which linearly dependent on the generated electricity and  No2 fluxes due to the wind in $f = \\{east, west, north, south \\}$ extremity of Puerto rucco, to compute an estimation of the variation of No2 total quantity due to these two factors:\n\n$  \\Delta C_{no2} = \\sum_{p \\in Power Plant}E_p(T)-(F_{east}+F_{west}+F_{north}+F_{south})*T $\n\nwhere $E_p(T)=P_p*T*EF_p$ the emission of the power plant $p$ over the periode of time $T=1$month, $F_{f}=\\sum_{i\\in f}C_i\\vec{U}_i\\cdot\\vec{e}_f$ are the average fluxes over the periode $T$. \n\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"NO2_concentration_month = np.zeros(12)\nestimated_NO2 = np.zeros(12)\ndate = []\n\n#  transform polor coordinate to (x,y)\nd = 1\nsigma = 0. # we do not need add the gaussian filter to remove noise, \n          # since we are averaging about 30 images for each periode\nfor i in range(0,12):\n    start_date = debut_date_list[i]\n    end_date = debut_date_list[i+1]\n    T_m = (end_date-start_date).total_seconds()\n    im_no2 = nan_filter_mean(folder_s5p_no2,images_s5p,start_date,end_date,band=2,sigma=sigma,acceptable_nan=0.1)\n\n    \n    im_no2 = im_no2\n    no2_north_face = np.mean(im_no2[-d:,:],axis=0)\n    no2_south_face = np.mean(im_no2[:d,:],axis=0)\n    no2_west_face = np.mean(im_no2[:,:d],axis=1)\n    no2_east_face = np.mean(im_no2[:,-d:],axis=1)\n\n    \n    mean_u = nan_filter_mean(folder_gfs,images_gfs,start_date,end_date,band=4,sigma=sigma)\n    mean_v = nan_filter_mean(folder_gfs,images_gfs,start_date,end_date,band=5,sigma=sigma)\n    mean_u = mean_u\n    mean_v = mean_v\n    speed_north_face = np.mean(mean_v[-d:,:],axis=0)\n    speed_south_face = np.mean(mean_v[:d,:],axis=0)\n    \n    speed_east_face = np.mean(mean_u[:,-d:],axis=1)\n    speed_west_face = np.mean(mean_u[:,:d],axis=1)\n    \n    flux_north = dx*(speed_north_face*no2_north_face).sum()\n    flux_south = -dx*(speed_south_face*no2_south_face).sum()\n    flux_east = dy*(speed_east_face*no2_east_face).sum()\n    flux_west = -dy*(speed_west_face*no2_west_face).sum()\n \n    month = start_date.ctime().split(' ')[1]\n    date.append(month)\n    estimated_NO2[i]  = thermal_power_plants['estimated_emission_mol_'+month].sum()-(flux_north+flux_south+flux_east+flux_west)*30*24*3600\n    \n    NO2_concentration_month[i] = dx*dy*im_no2.sum()\n    \n     \n    \n    \nNO2_concentration_month_normalize = NO2_concentration_month/NO2_concentration_month.mean()\nestimated_NO2_normalize = estimated_NO2/estimated_NO2.mean()\n\nplt.plot(date,NO2_concentration_month_normalize,date,estimated_NO2_normalize)\nplt.legend(['No2 concentration','estimated_NO2_normalize'])\nprint('the correlation between NO2 concentration and the generated electricity is ',np.corrcoef(NO2_concentration_month,estimated_NO2_normalize)[0,1])\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is hard to say that including the wind has a clear impact in correlation, since the correlation coefficent has \nonly increased from 0.42 to 0.5. But, we should keep in mind that the estimated generated electricity is not an accurate estimation. \nIn the next section, we will develop a physical model to predict the emission per power plant.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# 2. Modeling of NO2 concentration and emission factor\n\nNotation:\n\n- Np: number of power plant\n- i : power plant indice\n- $\\mathbf{x}_i$:  (x,y) position of the power plant in a 2D plane coordinate system\n- $ef_i$ :  emission factor of the power plant i (mol/Kwh)\n- $P_i$  : Power plant power (Mw)\n- $e_i$  : emision rate  $1./(3600)ef_p*1.e3*P_p$  (mol/s)\n- $c(t,\\mathbf{x})$ : No2 concetration at time t and position x ($mol/m^2$)\n- $c^{gfs}(t=d)$ : No2 concentration in day d from S5P data \n- $D^m$ :  molecular diffusivity coeffision of NO2  $m^2/s$\n- $\\Omega$ :  domain space\n- $\\delta \\Omega$: boundaries of the domain space\n\n## 2.1. NO2 modeling\nWe can model the variation of NO2 concentration by a convection-diffusion equation with source terms (all activities that can generate NO2 and reactions of NO2 with other chemical components). In the following, we will only consider power plant activities as the only source term for NO2:\n\n$$\\partial_t c+\\nabla\\cdot\\mathbf{u}c = D^m\\nabla^2{c}+\\sum_ie_i \\delta_{\\mathbf{x}_i},\\\\\n\\text{and initial solution:}\\,\\,\\,  c(t=t_0,x)=c(t_0,x)$$\n### 2.1.1 Boundary conditions\nSince we are working in a limited space, boundary conditions are needed to solve the above differential equation with inital solution. \nUnfortenately, we do not know the solution over all the time on the boundaries. The wind is mainly comming form the east and that were the in-flux of No2 will be important, it will be recommanded to do more analysis on No2 profile over y-axis in the east face in order to give a reliable boundary condition of No2 concentration along the East face. In the following, we will consider a simple estimation of dirichlet:\n\n$$c(t^{'},\\mathbf{x}) = c^{gfs}(t=d,\\mathbf{x})\\,\\, \\forall \\mathbf{x} \\in \\delta \\Omega\\,\\, \\text{and}\\,\\, d<t'<d+1day$$\n\n<u>Remark: The boundary conditions are the weakest part of this model, I will try to improve it in future work.</u>\n\n### 2.1.2 Turbulence and wind modeling\nThe model described below assumes an enough resolution with accurate speed values of the wind over the time. However, we only have a low resolution of wind map in periods of 6 hours. So we will need to overcome these two problem of space and time low resolution:\n1. Wind variation over time: the direct solution for that is to solve transport differential equation of the two component of the wind speed $\\mathbf{u}(t,\\mathbf{x})=(u,v)(t,\\mathbf{x})$. However, We are not going to this complicated solution. From the previous annalysis, we see that the wind does not vary it too much during the day. Therefore, we will use the wind data that we have for each 6 hours and interpolate to estimate the wind speed over all time. \n2. Low resolution: we have very rough resolution of the wind speed and also the initial solution, and even if we have a high resolution of these two parameters, it will be very CPU time consuming to solve numerically this differential equation in a large area with high resolution. We need therefore to consider a subgrid model, which will be an additional turbulent diffusion effect ([Turbulence Eddy Diffusivity](https://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-820-turbulence-in-the-ocean-and-atmosphere-spring-2007/lecture-notes/ch13.pdf)). The new equation can be writen as follows:\n\n    $$\\partial_t c+\\nabla\\cdot\\mathbf{u}c = D^m\\nabla^2{c}+\\nabla\\cdot(\\mathbf{\\varepsilon}\\cdot\\nabla(c))+\\sum_ie_i \\delta_{\\mathbf{x}_i}$$,\n    \n$\\mathbf{\\varepsilon}=(\\epsilon_{i,j})$ is the eddy diffusivity tensor. In the following, we consider an homogenouse and isotropic turbulence and we Prandlt mixing length model: $D^T=\\varepsilon_{x,x}=\\varepsilon_{y,y}\\approx L_m^2||\\mathbf{\\nabla{u}}||$, we take $L_m=C_smean(\\Delta x^{wind}, \\Delta y^{wind})$,where $C_s\\sim0.1-1$ is an hyperparameter to be tuned and $(\\Delta x^{wind}, \\Delta y^{wind})$ are the wind gfs real resolution. So we can simply write:\n\n$$\\partial_t c+\\nabla\\cdot\\mathbf{u}c = D\\nabla^2{c}+\\sum_ie_i \\delta_{\\mathbf{x}_i}$$,\n\n\nwhere $D=D^{m}+D^T$.\n\n$D^T\\sim 500m^2/s>>D^m$, so we can neglect the molecular diffusion.\n\n\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"\n# specify the time slot\nstart_time = datetime.datetime(2018, 9, 23,12,0,0) \nend_time = start_time+  datetime.timedelta(days=5)\n\n\n\ntime_s5p_list = get_time_stamp_list_s5p(s5p_filenames)\ntime_gfs_list = get_time_stamp_list_gfs(gfs_filenames)\n\n\ntimestamp_no2,im_no2 = get_data_in_time_slot(images_s5p[2],time_s5p_list,start_time,end_time)\n\n\n\ntimestamp_no2,im_no2 = remove_bad_images(timestamp_no2,im_no2,nan_limit_ratio=0.2)\n\n# We fix start_time as the time origine \n\nt_no2 = np.array([(t- start_time).total_seconds() for t in timestamp_no2])\n\nt_gfs = np.array([(t- start_time).total_seconds() for t in time_gfs_list])\n\n\n\n# gfs contains wind speed each 6hours and in 2D map of 'real resolution (4,9)'\n# continuous_wind_speed Object allows:- interpolate the (4,9) wind speed map to (148,475)\n#                                     - get wind speed-map in continuous time by interpolating over time\nwind = continuous_wind_speed(t_gfs,images_gfs[4],images_gfs[5],X,Y,(4,9))\n\n#  ---- downsample the resolution to by factor  reduction_factor\nreduction_factor = [4,5]\nnew_shape = (int(shape0[0]/reduction_factor[0]),int(shape0[1]/reduction_factor[1]))\nwind.downsample(new_shape)\nim_no2 = np.array(  [rebin(im, new_shape) for im in im_no2])\n    \n    \nX = np.pi/180*Rt*np.cos(lat*np.pi/180)*np.linspace(lon_limits[0],lon_limits[1],new_shape[1])\nY = np.pi/180*Rt*np.linspace(lat_limits[0],lat_limits[1],new_shape[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dx = np.abs(X[1]-X[0])\ndy = np.abs(Y[1]-Y[0])\nxy_loc, elec_year,emission_monthly = power_plant_location(thermal_power_plants,mean_im_no2_year.shape,lat_limits,lon_limits)\n\nXYp = [[x/reduction_factor[1],y/reduction_factor[0]]   for [x,y] in xy_loc]\nEm_p = [em['Nov']/(dx*dy)  for em in emission_monthly]\nm = 3\nfig,ax=  plt.subplots(m-1,3,figsize=(60,40))\nfor s in range(1,m):\n    # initilize the no2 solver\n    no2 = no2_diffusion_convection_solver(im_no2[s-1],X,Y,t_no2[s-1],t_no2[s],XYp, Em_p,D=0.6)\n    \n    bound = no2.boundary_condition(im_no2[s-1])\n    \n    no2.simulate_diffusion_convection_source_terms(wind,Cs=0.8)\n    \n    \n    #----------------- Plot the solution-------------#\n    im = ax[s-1,0].imshow(no2.C0[:,:],cmap='jet')\n    ax[s-1,0].set(title='initial concentration at t= '+str(timestamp_no2[s-1]))\n\n    cbar = fig.colorbar(im,ax=ax[s-1,0], format='%1.2e')\n    cbar.ax.tick_params(labelsize=25)\n    im =  ax[s-1,1].imshow(no2.C[:,:],cmap='jet')\n    ax[s-1,1].set(title='simulated  no2 concentration at t= '+str(timestamp_no2[s]))\n    cbar = fig.colorbar(im,ax=ax[s-1,1], format='%1.2e')\n    cbar.ax.tick_params(labelsize=25)\n    im = ax[s-1,2].imshow(im_no2[s],cmap='jet')\n    ax[s-1,2].set(title=' no2 concentration at t= '+str(timestamp_no2[s]))\n    cbar = fig.colorbar(im ,ax=ax[s-1,2], format='%1.2e')\n    cbar.ax.tick_params(labelsize=25)\n    #---------------------------------------------------#","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see the high concentration in north for both the simulated solution (using the estimated emission rate from the first part) and gfs image at $t=t_{end}$. Many patterns are common for the two images (center and rigth one), but it is still very hard to say that the two images matched prefectely. We will need to model the emission rate to get a closer prediction.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## 2.2 Emission factor modeling\n\n### 2.2.1 Decompose the solution to a combination of basis-space functions\nOur differential equation is linear, so we can apply the superposition principle. First, we decompose the\nsystem of Equation to two systems:\n$$ \\partial_t \\overline{c}+\\nabla\\cdot\\mathbf{u}\\overline{c} = D\\nabla^2{\\overline{c}},\\\\\n\\overline{c}(t=t_0,x)=c(t_0,x) \\\\\n\\text{and}\\, \\overline{c}(t^{'},\\mathbf{x}) = c^{gfs}(t=d,\\mathbf{x})\\,\\,\\, \\forall\\, \\mathbf{x} \\in \\delta \\Omega\\,\\, \\text{and}\\, d\\leq t'\\leq d+1day$$\n\n\nand \n\n$$ \\partial_t \\tilde{c}+\\nabla\\cdot\\mathbf{u}\\tilde{c} = D\\nabla^2{\\tilde{c}}+\\sum_ie_i \\delta_{\\mathbf{x}_i},\\\\\n\\tilde{c}(t=t_0,x)=0 \\\\\n\\text{and}\\, \\tilde{c}(t^{'},\\mathbf{x}) = 0\\,\\,\\, \\forall\\, \\mathbf{x} \\in \\delta \\Omega\\,, t\\geq t_0$$\n\nThe first system of equation can now be completely solved without the knoweledge of the emission rate of each power plant $(e_i)_{0,...,Np-1}$. However, the second system requires to have these emission rates. In fact, it is our main objectif to find these variables. The set of all possible solutions of the second system, when we vary the $(e_0,...,e_{Np-1})$ over $R^{Np}$, is a vectorial space of dimension $Np$. So, it will be very usefull to find a basis of this vectorial space. A trivial and a very usefull one can be found by solving for each $i\\in\\{0,...,Np-1\\}$ the following system:\n$$ \\partial_t \\tilde{c}_i+\\nabla\\cdot\\mathbf{u}\\tilde{c}_i = D\\nabla^2{\\tilde{c}_i}+e_i \\delta_{\\mathbf{x}_i},\\\\\n\\tilde{c}_i(t=t_0,x)=0 \\\\\n\\text{and}\\, \\tilde{c}_i(t^{'},\\mathbf{x}) = 0\\,\\,\\, \\forall\\, \\mathbf{x} \\in \\delta \\Omega\\,$$\n\n\nWe can decompose the final solution to n+1 contributions, $c=\\overline{c}+\\sum_ie_i\\tilde{c}_i$. \n\nSome remarks:\n1. only $\\overline{c}$ depends on the boundary conditions (remember it is our weak point in this model),\n2. $\\overline{c}$ does not depends on the emission,\n3. $\\tilde{c}_i$ for each i depends only on the location of the power plant i and we can expect that if two power plants are very close their respective solutions should be close to each other. We will also have some power plants that are in the same Finite Volume cell, so the numerical solution will be exactly the same in this case.\n\n\n### 2.2.2 Leat square fit to model the emission \nIn S5P data, we have the concentration for each day and with gfs.  We can numerically solve the system of equation of the unkown $C = (\\overline{c}, \\tilde{c}_0,...,\\tilde{c}_{Np-1})$ functions given above, with the initial solution $c^{gfs}(t=d)$ solution at day $d$  and compute an approximation solution at day $d+1$ $C(t=d+1)=(\\overline{c}(t=d+1), \\tilde{c}_0,...,\\tilde{c}_{Np-1}(t=d+1))$. Then, we can find the best parameter of the emission vairable $e_i$ by least square of:\n\n$$min(||c^{gfs}(t=d+1,x)-(\\overline{c}(t=d+1,x)+\\sum_ie_i\\tilde{c}_i(t=d+1))||_2),\\\\\n\\text{subject to}\\,\\,\\, e_i\\geq 0\\,\\,\\forall i\\in[0,Np-1]$$\n\n\n\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# we keep only power plant with important capacity,to avoid a kind of an overfitting of low power plant capacities that are close \n# to other sources of No2\npowerplant_filtered_capacity = thermal_power_plants[thermal_power_plants.capacity_mw>50.].reset_index()\nxy_p, el_p,em_p = power_plant_location(powerplant_filtered_capacity,shape0,lat_limits,lon_limits)\nname_p = np.array(powerplant_filtered_capacity['name'].tolist())\nxy_p = np.array(xy_p)\nel_p = np.array(el_p)\nem_p = np.array(em_p)\n# In case of more than one power plants are close to each ||xy_loc[i]-xy_loc[j]||_2<xy_dist_threshold, we take the barycenter of \n# these close power plant weigthed by their respective capacity\n\nxy_dist_threshold = 4\nXY_reduced = []\nEp_reduced = []\nname_reduce = []\ni = 0\nwhile(len(xy_p)>i):\n    delta_x = np.array(xy_p[i,:])-xy_p\n    I =  np.sum(delta_x**2,axis=1)<xy_dist_threshold\n    XY_reduced.append(np.mean(el_p[I]*xy_p[I,:],axis=0)/np.sum(el_p[I]))\n                      \n    Ep_reduced.append({m:  np.sum(np.array([e[m] for e in em_p[I]])) for m in months_list})\n    name_reduce.append(name_p[I])\n    I[i] = False  # we keep \n    xy_p = xy_p[~I]\n    em_p = em_p[~I]\n    el_p = el_p[~I]\n    name_p = name_p[~I]\n    i = i+1\nprint('The final number of source positions (power plant positions or barycenter of some power plant positions) that we are going to use: ',len(XY_reduced))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We keept only 9 power plant of capacity more than 50MW. After that, we keep only the barycenter of closer power plant. The final result is 8 positions to be using as location for emission in our simulation:\n\n- 7 position out 8 are real power plant positions\n- one position corresponds to the barycenter (weighted by the generate energy) of two real power plants","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Simulation of Basis functions ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"c_basis = []\n\nxy_p = [[x/reduction_factor[1],y/reduction_factor[0]]   for [x,y] in XY_reduced]  # reduce xy-index PP positions by the same factor as for no2 2d-map\n\n\nno2 = no2_diffusion_convection_solver(im_no2[0],X,Y,t_no2[0],t_no2[1],[], [],D=0.6)\n\nbound = no2.boundary_condition(im_no2[0])\n\nno2.simulate_diffusion_convection_source_terms(wind,Cs=0.9)\n\nc_bar = no2.C\n\nfor xyp in xy_p:\n    c0 = np.zeros(im_no2[0].shape)\n    no2 = no2_diffusion_convection_solver(c0,X,Y,t_no2[0],t_no2[1],[xyp], [1.],D=0.6)\n\n    bound = no2.boundary_condition(c0)\n\n    no2.simulate_diffusion_convection_source_terms(wind,Cs=0.9)\n    c_basis.append(no2.C)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure()\n\nplt.imshow(c_bar,cmap='jet')\nplt.title('$\\overline{c}$')\nplt.colorbar(format='%1.2e')\n#ax.set(title='$\\overline{c}$')\nfig,ax=  plt.subplots(len(c_basis)//2,2,figsize=(20,16))\nfor i in range(len(c_basis)):\n    ix = i%2\n    iy = i//2\n    im = ax[iy,ix].imshow(c_basis[i],cmap='jet')\n    ax[iy,ix].set(title='basis element: $c_'+str(i)+'$')\n    fig.colorbar(im, ax=ax[iy,ix],format='%1.2e')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Estimation of emission rate:**\n\n\n\nWe need find the best positive emission rate coefficient $(e_i)_{i=0..7}$ to minimize the erreur between the simulate solution $(\\overline{c}(t=d+1,x)+\\sum_ie_i\\tilde{c}_i(t=d+1))$ and our ground truth solution $c^{gfs}(t=d+1,x)$. In the current implementation the erreur is L2-norm of the delta between the two maps:\n\n$$min(||c^{gfs}(t=d+1,x)-(\\overline{c}(t=d+1,x)+\\sum_ie_i\\tilde{c}_i(t=d+1))||_2),\\\\\n\\text{subject to}\\,\\,\\, e_i\\geq 0\\,\\,\\forall i\\in[0,Np-1]$$\n\nRemark: after running some multiple simulations, I found that this fit is not working well all the time. The main problem was with $\\overline{c}$, as I pointed before $\\overline{c}$ depends on the boundary conditions and because we have a very modest model for the boundary conditions we can not trust its use on the fit, so I will add a coefficient $\\lambda$ to be included on the fit:\n\n$$min(||c^{gfs}(t=d+1,x)-\\big(\\lambda\\overline{c}(t=d+1,x)+\\sum_ie_i\\tilde{c}_i(t=d+1))\\big)||_2),\\\\\n\\text{subject to}\\,\\,\\, e_i\\geq 0\\,\\,\\forall i\\in[0,Np-1]\\,\\,\\, \\text{and}\\,\\,\\,\\lambda\\geq \\lambda_{min}$$\n\n\n<u>In the following, we take $\\lambda_{min}=0.5$. This is an ad-hoc value to overcome this problem, we will need to improve the boundary conditions in future.</u>\n","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Since, we neglect lot of other source term of No2 from other different activities, it will be very wise to focus the least square only on area near to the simulate power plant locations. We will then use a mask to keep only the region of interest:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sh = c_basis[0].shape\nXint,Yint = np.meshgrid(np.linspace(0,sh[1]-1,sh[1]),np.linspace(0,sh[0]-1,sh[0]))\ndr = np.array([np.sqrt((Xint-xy[0])**2+ (Yint-xy[1])**2) for xy in XYp  ])\nroi = np.min(dr,axis=0)<7.\nplt.imshow(roi)\nplt.title('Region of Interest')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"coef,no2_fitted = emission_factor_lsq_with_lambda_nnls(c_bar,c_basis,im_no2[1],roi,lambda_lim=0.5)\n\nfig,ax=  plt.subplots(1,2,figsize=(40,14))\nim = ax[0].imshow(no2_fitted,cmap='jet',vmin=0,vmax=2.5e-5)\nax[0].set(title='fitted solution')\ncbar =fig.colorbar(im,ax=ax[0],format='%1.2e')\ncbar.ax.tick_params(labelsize=25)\n\nim = ax[1].imshow(im_no2[1],cmap='jet',vmin=0,vmax=2.5e-5)\nax[1].set(title='ground true soltution')\ncbar = fig.colorbar(im,ax=ax[1],format='%1.2e')\ncbar.ax.tick_params(labelsize=25)\nprint('emission rate coefficient',coef)\n\n\n#  delta between the fitted solution and the reference solution\n\nfig,ax = plt.subplots(1,2,figsize=(20,10))\nim = ax[0].imshow(np.abs(no2_fitted-im_no2[1])*roi,cmap='jet',vmax=3.e-5)\nax[0].set(title='delta between the fitted solution and groundtruth on region of interest')\nfig.colorbar(im,ax=ax[0],format='%1.2e')\nim = ax[1].imshow(np.abs(no2_fitted-im_no2[1]),cmap='jet',vmax=3.e-5)\nax[1].set(title='delta between the fitted solution and groundtruth')\nfig.colorbar(im,ax=ax[1],format='%1.2e')\n\ndelta = np.abs(no2_fitted-im_no2[1])\nprint(' relative std error in ROI ',np.std(delta[roi])/mean_im_no2_year.mean())\nprint(' relative std error in all the domain ',np.std(delta)/mean_im_no2_year.mean())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now let's run this model for different random days over the year:\n - First, we remove bad no2 data: (no2 images with more than 20% of nans)\n - We simulate the concentration between two successive good available S5P data with condition $\\Delta T <1.5 day$\n - Each simulation consists in simulating $\\overline{c}$ and the 8 basis function $\\tilde{c}_{i=0..,7}$\n ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"columns = ['date','std_erro','std_erro_roi'] + ['emission_rate_'+str(i)  for i in range(8)]\n\nprint(columns)\npd_sim = pd.DataFrame(columns=columns)\n\nstart_time = datetime.datetime(2018, 7, 8,12,0,0)  \nend_time = start_time+  datetime.timedelta(days=365)\n\n\n\ntime_s5p_list = get_time_stamp_list_s5p(s5p_filenames)\ntime_gfs_list = get_time_stamp_list_gfs(gfs_filenames)\n\n\ntimestamp_no2,im_no2 = get_data_in_time_slot(images_s5p[2],time_s5p_list,start_time,end_time)\n\n\n\ntimestamp_no2,im_no2 = remove_bad_images(timestamp_no2,im_no2,nan_limit_ratio=0.2)\n\n# We fix start_time as the time origine \n\nt_no2 = np.array([(t- start_time).total_seconds() for t in timestamp_no2])\n\nt_gfs = np.array([(t- start_time).total_seconds() for t in time_gfs_list])\n\n\n\n# gfs contains wind speed each 6hours and in 2D map of 'real resolution (4,9)'\n# continuous_wind_speed Object allows:- interpolate the (4,9) wind speed map to (148,475)\n#                                     - get wind speed-map in continuous time by interpolating over time\nwind = continuous_wind_speed(t_gfs,images_gfs[4],images_gfs[5],X,Y,(4,9))\n\n#  ---- downsample the resolution by factor of reduction_factor to speedup the simulation\nreduction_factor = [4,5]\nnew_shape = (int(shape0[0]/reduction_factor[0]),int(shape0[1]/reduction_factor[1]))\nwind.downsample(new_shape)\nim_no2 = np.array(  [rebin(im, new_shape) for im in im_no2])\n\n\nX = np.pi/180*Rt*np.cos(lat*np.pi/180)*np.linspace(lon_limits[0],lon_limits[1],new_shape[1])\nY = np.pi/180*Rt*np.linspace(lat_limits[0],lat_limits[1],new_shape[0])\n\n\n\nxy_p = [[x/reduction_factor[1],y/reduction_factor[0]]   for [x,y] in XY_reduced]  # reduce xy-index PP positions by the same factor as for no2 2d-map\n\n\nim_no2[0][im_no2[0]<=1.e-8] = 1.e-8  # remove negative values\nfor s in range(1,len(im_no2)):\n\n    if((t_no2[s]-t_no2[s-1])<3600*24*1.5):\n        no2 = no2_diffusion_convection_solver(im_no2[s-1],X,Y,t_no2[s-1],t_no2[s],[], [],D=0.6)\n\n        bound = no2.boundary_condition(im_no2[s-1])\n\n        no2.simulate_diffusion_convection_source_terms(wind,Cs=0.9)\n\n        c_bar = no2.C\n        c_basis = []\n        for xyp in xy_p:\n            c0 = np.zeros(im_no2[s-1].shape)\n            no2 = no2_diffusion_convection_solver(c0,X,Y,t_no2[s-1],t_no2[s],[xyp], [1.],D=0.6)\n\n            bound = no2.boundary_condition(c0)\n\n            no2.simulate_diffusion_convection_source_terms(wind,Cs=0.9)\n            c_basis.append(no2.C)\n        im_no2[s][im_no2[s]<=1.e-8] = 1.e-8   # remove negative values\n        emission_coeff,no2_fitted = emission_factor_lsq_with_lambda_nnls(c_bar,c_basis,im_no2[s],roi,lambda_lim=0.5)\n\n        delta = np.abs(no2_fitted-im_no2[s])\n        row_pd = {'emission_rate_'+str(i):emission_coeff[i]*(no2.dx*no2.dy)/(t_no2[s]-t_no2[s-1]) for i in range(8)}\n        row_pd['date'] = timestamp_no2[s-1]\n\n        row_pd['std_erro'] = np.std(delta)/np.mean(im_no2[s])\n        row_pd['std_erro_roi'] = np.std(delta[roi])/np.mean(im_no2[s][roi])\n\n        pd_sim = pd_sim.append(row_pd, ignore_index=True)\n        print(' ')\n        print(' relative std error in ROI ',np.std(delta[roi])/mean_im_no2_year.mean())\n        print(' relative std error in all the domain ',np.std(delta)/mean_im_no2_year.mean())\n\n\n\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Daily emission rate per power plant","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(10,8))\npd_sim_filtered = pd_sim[pd_sim.std_erro_roi<0.5]\n\npd_sim_filtered.plot(x='date',y=['emission_rate_'+str(i) for i in range(8)], ax=ax)\nax.legend(['emission rate of '+str(name_reduce[i])+'power plant' for i in range(8)])\nax.set(ylabel='mol/s')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Monthly emission rate per power plant","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pd_sim_filtered['month'] = pd_sim_filtered.date.apply(lambda r:r.month)\npd_sim_monthly = pd_sim_filtered.groupby('month').mean().reset_index()\npd_sim_monthly","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(10,8))\n\npd_sim_monthly.plot(x='month',y=['emission_rate_'+str(i) for i in range(8)], ax=ax,marker='.')\n\nax.legend(['emission rate of '+str(name_reduce[i])+' power plant' for i in range(8)])\nax.set(ylabel='mol/s')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3- Parameter study and validation\nIn this final part, we study the sensitvity of the model to the different hyperparmeters and the validation of the moedl:\n\n\n**to be continued**","execution_count":null}],"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}