{"cells":[{"metadata":{},"cell_type":"markdown","source":"#
N02 PANACHE FITTING
\n
Louis Bouvier & Pablo Richard
\n

\n
Kaggle competition: DS4G - Environmental Insights Explorer
\n

\n
Exploring alternatives for emissions factor calculations
"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false,"scrolled":true},"cell_type":"code","source":"%%javascript\n// little fix to enumerate equations in markdown\n// from https://github.com/ipython-contrib/jupyter_contrib_nbextensions/tree/master/src/jupyter_contrib_nbextensions/nbextensions/equation-numbering\n\nMathJax.Hub.Config({\n TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n});","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nfrom matplotlib import ticker\nimport folium\nimport os\nimport rasterio as rio\nimport numpy as np\nimport pandas as pd\nimport tifffile as tiff \nfrom IPython.display import Image\nfrom datetime import datetime\nimport seaborn as sns\nimport matplotlib.cm as cm \nfrom folium import plugins \nfrom datetime import datetime\nfrom matplotlib import ticker\nfrom scipy import stats\nfrom scipy.special import erf","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## I - Introduction and overview of the notebook"},{"metadata":{},"cell_type":"markdown","source":"### Goal\n\nDevelop a methodology to calculate an average historical emissions factor of electricity generated for a sub-national region, using remote sensing data and techniques.\n(from competition overview)\n\n\n### Emission factor\n\nAn emissions factor is a representative value that attempts to relate the quantity of a pollutant released to the atmosphere with an activity associated with the release of that pollutant. These factors are usually expressed as the weight of pollutant divided by a unit weight, volume, distance, or duration of the activity emitting the pollutant (*e.g.*, kilograms of particulate emitted per megagram of coal burned). Such factors facilitate estimation of emissions from various sources of air pollution. In most cases, these factors are simply averages of all available data of acceptable quality, and are generally assumed to be representative of long-term averages for all facilities in the source category (*i.e.*, a population average).\n\nThe general equation for emissions estimation is:\n\nE = A x EF x (1-ER/100)\n\nwhere:\n\n E = emissions;\n A = activity rate;\n EF = emission factor, and\n ER = overall emission reduction efficiency, %\n\nfrom [[1]](#References)"},{"metadata":{},"cell_type":"markdown","source":"Since with satellite data we have a way to estimate emissions, and the plant dataset highlights some electricity generation assessments, we can deem the emission factor in the following way :\n\n$$EF \\approx E/A$$"},{"metadata":{},"cell_type":"markdown","source":"#### Let's have in mind some rough estimates of NOx emissions per primary fuel"},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Image(filename='/kaggle/input/external-data/emission_per_pf.png')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"source : [[2]](#References)"},{"metadata":{},"cell_type":"markdown","source":"Since we have access to NO2 satellite data - which can be linked to NOx emissions thanks to mass conversion regarding a chemical transformation according to [[3]](#References) - we are going to focus on some specific primary fuels to deem power plants emissions. \n\nTo do so, one can notice that the main sources (among those used in Puerto Rico, see below please) of emissions per MWh of electricity produced are coal, gas, oil and solar. We will thus exclude hydropower and wind for a first approximation. Indeed the 10-100 factor between those two and the rest seems enough to discard them. This choice will be bolstered by the comparison of electricity generation below : coal, gas and oil account for the majority sources of electricity."},{"metadata":{},"cell_type":"markdown","source":"##
Content of the rest of this notebook
"},{"metadata":{},"cell_type":"markdown","source":"## [II.](#II---Data-exploration) - Data exploration\n\n ### A) A glimpse at initial data\n \n ### B) Modification of the *GPPD* and justification\n #### a) Link between capacity and energy generation \n #### b) Add information from another dataset specific to US including Puerto Rico \n \n \n## [III.](#III---First-attempt-to-compute-an-emission-factor:-focus-on-clusters-and-apply-gaussian-filters) - First attempt to compute an emission factor: focus on clusters and apply gaussian filters\n\n ### A) NO2 data pre-processing\n #### a) Loading and plotting\n #### b) Selection with *NaN* threshold\n #### c) Data interpolation - 2D linear per frame using *pandas.interpolate*\n \n ### B) Start with one cluster\n \n ### C) From instantaneous NO2 to plant emission factor\n \n ### D) Extend approach to clusters and sum to get a first estimation for the whole island\n\n\n## [IV.](#IV---Second-attempt:-prepare-a-broader-dataset-with-horizontal-wind,-temperature-and-relative-humidity) - Second attempt: prepare a broader dataset with horizontal wind, temperature and relative humidity\n\n ### A) Extract *GFS* 2D fields and match the dates with NO2\n \n ### B) From 2D fields to 8 1D fields\n\n\n## [V.](#V---Fitting:-physical-analysis-to-derive-an-analytical-1D-model) - Fitting: physical analysis to derive an analytical 1D model\n\n ### A) Physical discussion for the mathematical form of the surfacic NO2 density\n \n ### B) Let's go into an analytical 1D model\n \n ### C) Fitting implementation\n\n\n## [VI.](#VI---Final-results-and-critics) - Final results and critics\n\n ### A) Results and answer to the competition\n \n ### B) Comments and perspectives\n\n\n## [References.](#References)"},{"metadata":{},"cell_type":"markdown","source":"## General ideas of this submission:\n\nTo estimate the year averaged emissions factor we first preprocess data, modifying the power plants database with physics arguments and adding monthly information for various power-plants with an external dataset, and interpolating NO2 pictures in a certain limit of number of *NaN* values. This step is crucial to be able to estimate emission factors at a smaller slice than by year.\n\nThen, a first \"brutal\" approach is taken: we group NO2-emitting power-plants in clusters regarding their locations (since some of them are close to each others) and we define a certain mask to locally extract NO2 above those clusters. At this step, we have merely mixed *GPPD* and *S5P_NO2* datasets to provide an average historical emissions factor of electricity generated for a sub-national region, using remote sensing data and techniques. We discuss the limits of this first approach in the end of section III.\n\nTo improve those results, we refine our approach: we merge NO2 pictures with the *GFS* dataset to get wind, temperature and humidity information. In addition, motivated by a physics background, we create a new dataset which presents information in a different way: instead of working with 2D fields, we extract 1D density fields along discrete directions and in a certain perimeter for each cluster of plants and each accessible day. For each cluster, we then have per day 8 1D fields (one for each main discrete cardinal point) for temperature, humidity, wind and NO2. \n\nWe then define an additve model for NO2 sources and diffusion and fit in the discrete direction that corresponds to the wind local direction (thus per cluster of power plants). Thanks to this model, we distinguish plant-related and background-related (nature, maybe cities and airports?) emissions per cluster of plants and **per day**. We believe this sharper focus is one way to deduce plant emissions from the total NO2 tropospheric field.\n\nWe then display final results on emission factors and comments. Having no reliable data for energetic generation over smaller time slices (even though trying to do it in first part), we were not able to derive instantaneous emission factors, neither marginal. However, we are fairly positive that adding these energetical production data to our daily emission estimation heuristic would be quite interesting!"},{"metadata":{},"cell_type":"markdown","source":"## II - Data exploration"},{"metadata":{},"cell_type":"markdown","source":"*Go back to table of contents [⤴](#Content-of-the-rest-of-this-notebook)*"},{"metadata":{},"cell_type":"markdown","source":"### A) A glimpse at initial data"},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"### functions created by the Kaggle community\n# Code source: https://www.kaggle.com/paultimothymooney/overview-of-the-eie-analytics-challenge \ndef plot_points_on_map(dataframe, begin_index, end_index,\n latitude_column, latitude_value, \n 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 str(df.latitude[i:i+1])+str(df.longitude[i:i+1]))\n folium.Marker([df[latitude_column].iloc[i],\n df[longitude_column].iloc[i]],popup=popup).add_to(plot)\n return plot\n\ndef overlay_image_on_puerto_rico(file_name,band_layer):\n band = rio.open(file_name).read(band_layer)\n m = folium.Map([lat, lon], zoom_start=8)\n folium.raster_layers.ImageOverlay(image=band,\n bounds=[[18.6,-67.3,],[17.9,-65.2]],\n colormap=lambda x: (1, 0, 0, x)).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,\n 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":{},"cell_type":"markdown","source":"#### Overview of the plants"},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"power_plants = pd.read_csv('/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\npower_plants = split_column_into_new_columns(power_plants,'.geo','latitude',50,66)\npower_plants = split_column_into_new_columns(power_plants,'.geo','longitude',31,48)\npower_plants['latitude'] = power_plants['latitude'].astype(float)\na = np.array(power_plants['latitude'].values.tolist()) # 18 instead of 8\npower_plants['latitude'] = np.where(a < 10, a+10, a).tolist() \n\nlat=18.200178; lon=-66.664513\nplot_points_on_map(power_plants,0,425,'latitude',lat,'longitude',lon,9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Details on power plants"},{"metadata":{"trusted":true,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Capacity = power_plants[\"capacity_mw\"]\nElectricity_generation = power_plants[\"estimated_generation_gwh\"]\nPrimary_fuels = power_plants[\"primary_fuel\"]\nfor p_f in Primary_fuels.unique():\n E_g = Electricity_generation[Primary_fuels == p_f].sum()\n C_g = Capacity[Primary_fuels == p_f].sum()\n print(\"Capacity = {:.2f} MW, Electricity generation in GWh = {:.2f} for primary fuel : {}\".format(C_g,E_g,p_f))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since electricity generation from solar primary fuel is more than 1000 times lower than from oil, coal and gas, we first focus only on the following primary fuels (for emissions estimations): oil, coal and gaz."},{"metadata":{},"cell_type":"markdown","source":"With the example of oil, can we emphasize a link between the estimated generation of electricity and the capacity of each plant?"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Oil_capa = np.array(Capacity[Primary_fuels == 'Oil'])\nOil_gen = np.array(Electricity_generation[Primary_fuels == 'Oil'])\nprint('Oil capacities (MW): \\n', Oil_capa)\nprint('Oil electricity generations (GWh): \\n', np.array(Oil_gen, dtype=np.float16))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We notice that the estimated electricity generation for those oil plants can be computed from the capacity, mutiplied by a constant factor (time unit).\n\nWe will be able to extend to smaller time slices with this assumption."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"# We load the time per primary fuel to get a lower bound of time \n# generation estimation with assumption of maximal capacity used\n\nTarget_primary_fuels = ['Oil','Gas','Coal']\n\nfor p_f in Target_primary_fuels:\n time_factor = np.array(1000*Electricity_generation[Primary_fuels == p_f])/np.array(Capacity[Primary_fuels == p_f])\n print('Time factor in hours for {}:\\n{}'.format(p_f, np.array(time_factor, dtype=np.float32)))\n\nprint('\\nNumber of hours in a year : ',365*24)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As one might see, in the last case of coal, the time needed to generates its anual energy using a maximal capicity exceeds a year!"},{"metadata":{},"cell_type":"markdown","source":"#### Let's go into details for the coal plant: is data consistent?"},{"metadata":{},"cell_type":"markdown","source":"It seems there is an issue of consistency between capacity and estimated electricity generation. Regarding data from other sources [[4]](#References), we can rely on capacity of the previous dataset more than on estimation.\n\nThis link between capacity and generation is described in [[5]](#References)."},{"metadata":{},"cell_type":"markdown","source":"### B) Modification of the *GPPD* and justification"},{"metadata":{},"cell_type":"markdown","source":"#### a) Link between capacity and energy generation "},{"metadata":{},"cell_type":"markdown","source":"Let's start with a physics assessment on thermodynamic systems and heat engines:\n\n$$ \\forall i,T \\in (I \t\\times \\mathbb R), \\,\\, Cmax_i*T = \\int_{0}^{T}Cmax_idt\\geq \\int_{0}^{T} C_i(t)dt = G_{T,i} $$ "},{"metadata":{},"cell_type":"markdown","source":"Where:\n- $I$ is a set of indexes of plants : each index corresponds to a plant. \n- $T$ is a time duration : a period during which we assess the production of a plant.\n- $Cmax_i$ is the maximum capacity (a constant) of the plant indexed by i (capacity_MW in our dataset) (measure : power).\n- $C_i(t)$ is the instantaneous capacity of the plant indexed by i at time t. It encapsulates to which extent the plant is running to produce electricity. If it is off, the value is 0 W. (measure : power).\n- $G_{T,i}$ is the estimated generation of the plant indexed by i during the time duration T (estimated_generation_GWh in our dataset) (measure : energy)."},{"metadata":{},"cell_type":"markdown","source":"*NB*: The time origin is translated to zero in the formula above to reduce writing, but the wear on plants and conditions may influence the instantaneous capacity."},{"metadata":{},"cell_type":"markdown","source":"#### Consequences on our dataset"},{"metadata":{},"cell_type":"markdown","source":"The idea here is to check whether the inequality above is respected for every plant in our Global Power Plant database :\n\ncapacity_mw.T >= estimated_generation_gwh*1000, where T is one year in our case."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"max_energy = power_plants[\"capacity_mw\"]*24*365\nestimated_energy = power_plants[\"estimated_generation_gwh\"]*1000","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"plt.bar(power_plants[\"primary_fuel\"],estimated_energy/max_energy)\nplt.title(\"Ratio estimation/max energy per primary fuel\")\nplt.xlabel(\"primary fuel\")\nplt.ylabel(\"ratio estimation/max energy\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As you can see some plants have a ratio greater than one, which means data is not consistent.\nWe are going to fix values for those having this ratio greater than one.\n\n- According to [[6]](#References) we must set a capacity factor (the estimation/max ratio) of $\\textbf{0.54}$ for coal plants in 2018 in the US.\n- According to [[7]](#References) we must set a capacity factor (the estimation/max ratio) of $\\textbf{0.32}$ for hydro plants."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"power_plants.loc[power_plants.primary_fuel == \"Coal\", \"estimated_generation_gwh\"] = \\\n0.54*power_plants[\"capacity_mw\"][power_plants[\"primary_fuel\"] == \"Coal\"]*24*365/1000\npower_plants.loc[power_plants.primary_fuel == \"Hydro\", \"estimated_generation_gwh\"] = \\\n0.32*power_plants[\"capacity_mw\"][power_plants[\"primary_fuel\"] == \"Hydro\"]*24*365/1000","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check again:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"max_energy = power_plants[\"capacity_mw\"]*24*365\nestimated_energy = power_plants[\"estimated_generation_gwh\"]*1000","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"plt.bar(power_plants[\"primary_fuel\"],estimated_energy/max_energy)\nplt.title(\"Ratio estimation/max energy per primary fuel\")\nplt.xlabel(\"primary fuel\")\nplt.ylabel(\"ratio estimation/max energy\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"*NB*: One can notice that the ratios for Gas and Oil are greater than in literature. Nonetheless we will keep them as such for the moment since we would need a dataset specific to Puerto Rico plants to adjust precisely. The goal here was first to avoid highly inconsistent data. We thus export."},{"metadata":{},"cell_type":"markdown","source":"#### b) Add information from another dataset specific to US including Puerto Rico "},{"metadata":{},"cell_type":"markdown","source":"Inspired by the valuable notebook [[8]](#References), we are going to complete our power plant dataset with month slices and information on energy generation and emissions."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"add_data = pd.read_csv('/kaggle/input/external-data/gppd_120_pr_ef.csv') #additional dataset","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"completed_gppd = pd.merge(power_plants, add_data, how='left', on=\".geo\") #merging with kaggle dataset\n\n# deleting '_x' at the end of labels\nfor column in completed_gppd.columns:\n if column[-2:] == '_x':\n completed_gppd.rename(index=str, columns={column: column[:-2]}, inplace=True)\n \n#we store our corrected dataset\ncompleted_gppd.to_csv(r'Power_plants_corrected.csv', index = False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We now load data from this corrected database for the rest of this processing notebook: "},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"power_plants = completed_gppd.copy()\nprint('Power plants dataset correction is done.')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## III - First attempt to compute an emission factor: focus on clusters and apply gaussian filters "},{"metadata":{},"cell_type":"markdown","source":"*Go back to table of contents [⤴](#Content-of-the-rest-of-this-notebook)*"},{"metadata":{},"cell_type":"markdown","source":"### A) NO2 data pre-processing"},{"metadata":{},"cell_type":"markdown","source":"We will need some spatial information on those satellite images to deem emissions linked to local sources: oil, gas and coal plants. We thus create a function to get the coordinates (longitude, latitude) of each pixel of an image thanks to its affine transformation."},{"metadata":{},"cell_type":"markdown","source":"*NB*: Since the spatial location we consider (the Puerto Rico island) is quite \"small\" compared to characteristic distances that require 3D computations, we will approximate the (longitude,latitude) coordinates as its projection in a 2D field over the island. \nThat is to say, we will neglect the influence of the ellipsoid shape of the earth on coordinates and distances computations below. This assumption does not seem absurd knowing pixel resolution of our datasets."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def coordinates(file_name):\n '''get all the coordinates of the pixels in the file'''\n data = rio.open(file_name)\n transform_coef = data.transform\n height, length = data.read().shape[1:]\n coord = np.array([[transform_coef*(i, j) for i in range(length)] for j in range(height)])\n return(coord)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### a) Loading and plotting"},{"metadata":{},"cell_type":"markdown","source":"In the following cells, we extract the tropospheric NO2 density field. This is motivated by the following ideas: \n\n\n- since NO2 has a short lifetime (see figures below) and troposphere is about 15 km high, molecules emitted at the level of the ground mostly won't have time to reach an altitude above.\n\n\n\n- according to papers like [[9]](#References), based on the same dataset from TROPOMI and according TROPOMI's documentation itself (see \"TROPOMI ATBD of the total and tropospheric NO2 data products\")[[10]](#References), the way NO2 is estimated entails more difficulties to distinguish human sources and natural sources in higher atmosphere (pressure and wind variations play no mean part in stratosphere).\n\n\n- According to TROPOMI documentation [[10]](#References) each image shows an instantaneous field of NO2 at a certain UTC time. We store this time to merge with relevant information in the *GFS* data. We will also need to extrapolate to deduce daily emissions from a single picture of instantaneous emission per day. See below please."},{"metadata":{},"cell_type":"markdown","source":"After a first attempt to extract NO2 data we faced two issues:\n- Some duplications in dates\n- Some *NaN* values in spatial data\n\nWe thus chose to interpolate pictures with a reasonable (thus small) proportion of *NaN* values and to deal with dates through the following cells."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def plot_scaled(file_name, title=None):\n '''plot a scaled colormap (avoid NAN issues)'''\n fig = plt.figure(figsize=(6,2))\n plt.title(title) if title is not None else None\n Y,X = np.meshgrid(np.arange(475),np.arange(148))\n vmin, vmax = np.nanpercentile(file_name, (5,95)) # 5-95% stretch\n img_plt = plt.pcolormesh(Y,X,file_name, vmin=vmin, vmax=vmax)\n plt.axis('equal')\n plt.colorbar()\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/'\n\nCoords_per_time = []#we store coordinates computed thanks to the previous function\ntropospheric_NO2_per_time_raw = []#we store the 2D scalar field\nTime_laps = []#we store dates \nHours_NO2_pictures = []#we store hour when satellite takes instantaneous picture\ndoublons_list = []\nextract_iteration = 0\nfor filename in os.listdir(path):\n extract_iteration += 1\n try:\n path_ = path+filename\n tropospheric_NO2_per_time_raw.append(rio.open(path_).read()[1])\n #print('Extracting file number {} named : {}'.format(extract_iteration, filename))\n date_begin = filename[8:16]\n date_end = filename[24:32]\n hour_picture = filename[17:23]\n if date_begin in list(map(lambda couple : couple[0], Time_laps)):\n #print('duplicated!')\n doublons_list.append(date_begin)\n Time_laps.append([date_begin,date_end])\n Hours_NO2_pictures.append(hour_picture)\n Coords_per_time.append(coordinates(path_))\n except:\n print('uncaught file : {}'.format(filename))\n \n\nTime_laps = np.array(Time_laps)\nCoords_per_time = np.array(Coords_per_time, dtype=np.float32)\ntropospheric_NO2_per_time_raw = np.array(tropospheric_NO2_per_time_raw, dtype=np.float16)\nHours_NO2_pictures = np.array(Hours_NO2_pictures)\n \n\nn_t, n_y, n_x = np.shape(tropospheric_NO2_per_time_raw)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see sometimes that two (never more) pictures have the same date, we call this a *doublon*. Since there is always one way better than the other, we just keep the best:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"print('Deleting duplicates')\nindex_to_del = []\nfor date in doublons_list:\n frame_list = []\n indexes = []\n for t in range(n_t):\n if Time_laps[t][0]==date:\n indexes.append(t)\n frame_list.append(tropospheric_NO2_per_time_raw[t])\n nan_list = [len(frame[np.isnan(frame)]) for frame in frame_list]\n for i, ind in enumerate(indexes):\n if nan_list[i]==max(nan_list):#when we see a duplicated date, we keep the one with less NaN values\n index_to_del.append(ind)\n\nindex_to_keep = [ind for ind in range(n_t) if not ind in index_to_del]\nn_t = len(index_to_keep)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#keep data after duplicates exclusion\nTime_laps = Time_laps[index_to_keep]\nCoords_per_time = Coords_per_time[index_to_keep]\ntropospheric_NO2_per_time_raw = tropospheric_NO2_per_time_raw[index_to_keep]\nHours_NO2_pictures = Hours_NO2_pictures[index_to_keep]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We create a physical timeline named *timestamp* using the Unix timestamp unity:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"timestamp = np.array([datetime.timestamp(datetime.strptime(str(int(date_couple[0])),\n '%Y%m%d'))\n for date_couple in Time_laps])\ntimestamp_hours = np.array([datetime.timestamp(datetime.strptime(str(int(date_couple[0]))+\\\n hour, '%Y%m%d%H%M%S'))\n for date_couple, hour in zip(Time_laps, Hours_NO2_pictures)])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We sort the data according to time:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"print('sorting data ...')\ntotal_data = {instant: (lap, hour, coords, frames)\n for instant, lap, hour, coords, frames\n in zip(timestamp,\n Time_laps,\n Hours_NO2_pictures,\n Coords_per_time,\n tropospheric_NO2_per_time_raw)}\n\ntimestamp = np.array(sorted(timestamp))\ntimestamp_hours = np.array(sorted(timestamp_hours))\nTime_laps = np.array([total_data[instant][0] for instant in timestamp])\nHours_NO2_pictures = np.array([total_data[instant][1] for instant in timestamp])\nCoords_per_time = np.array([total_data[instant][2] for instant in timestamp])\ntropospheric_NO2_per_time_raw = np.array([total_data[instant][3] for instant in timestamp])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Printing first 10 pictures, we see a lot of empty pixels due to *NaN* values in our frames:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"for t in range(10):\n plot_scaled(tropospheric_NO2_per_time_raw[t],\n title='Initial frame '+\\\n format(datetime.fromtimestamp(int(timestamp_hours[t])).strftime('%Y/%m/%d %H:%M:%S')))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### b) Selection with *NaN* threshold"},{"metadata":{},"cell_type":"markdown","source":"The following cells aim at interpolating NO2 pictures in a certain extent. Indeed, we will proceed to 2D interpolations on pictures that do not have too many *NaN* values (it would make no sense to infer values on a grid that has too many missing meaningful values). We thus first deem a threshold of amount of *NaN* values under which we interpolate. Pictures that have too many *NaN* values are simply excluded."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"information_trust = np.array(list(map(lambda frame : len(frame[np.isnan(frame)]),\n tropospheric_NO2_per_time_raw)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"plt.scatter(timestamp_hours, 100*information_trust/(n_x*n_y), s=2)\nplt.xlabel(r'Unix timestamp', fontsize=12)\nplt.ylabel(r'Percentage of unknown pixels $(NaN)$ - $(\\%)$', fontsize=12)\nplt.title(r'Noise percentage on frames against time', fontsize=18, pad=10)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"thresholds = np.array(list(range(np.max(information_trust))))\namounts = np.array([len(information_trust[information_trust<=threshold])\n for threshold in thresholds])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#calibrated by hand : maximum number of NaN accepted pixels in one picture (5%)\noptim_threshold = 2240\n\nplt.scatter(100*thresholds/(n_x*n_y),\n amounts,\n s=0.1,\n label='Cumulated histogram')\nplt.axvline(x=100*optim_threshold/(n_x*n_y),\n linestyle='--',\n linewidth=1,\n c='orange',\n label='Optimal threshold')\nplt.axhline(y=n_t, linestyle='--',\n linewidth=1,\n c='purple',\n label='Total number of pictures')\nplt.xlabel(r'Unknown pixels $(NaN)$ - $(\\%)$', fontsize=12)\nplt.ylabel(r'Number of frames with a NaN percentage smaller', fontsize=12)\nplt.legend(fontsize=15)\nplt.title(r'Finding the optimal $NaN$ threshold', fontsize=18, pad=15)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Hence in the following, we will be trying to interpolate the 2D frames with a number of *NaN* under this threshold."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"timestamp = np.array([timestamp[t] for t in range(n_t)\n if information_trust[t]2**, which is a great increase compared to the 129 without *NaN* values!\nLet's process this data."},{"metadata":{},"cell_type":"markdown","source":"### B) Start with one cluster"},{"metadata":{"scrolled":false,"trusted":false,"_kg_hide-input":true},"cell_type":"code","source":"### collect locations of gas, oil, and coal\nTarget_plant_coordinates = power_plants[[\"longitude\",\"latitude\"]][np.isin(power_plants[\"primary_fuel\"],Target_primary_fuels)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's see where the plants with fossil primary fuels are located "},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"lat=18.200178; lon=-66.664513\ntarget_power_plants = power_plants[np.isin(power_plants[\"primary_fuel\"],Target_primary_fuels)].copy()\nplot_points_on_map(target_power_plants,0,425,'latitude',lat,'longitude',lon,9)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's estimate a characteristic distance for the focus on a cluster and apply a conversion for longitudes/latitudes"},{"metadata":{},"cell_type":"markdown","source":"We are going to set a characteristic distance for gaussian filtering based on this assumption : NO2 molecules have a lifetime $\\tau$ that varies from 1.8 to 7.5 hours, regarding literature [[11]](#References) and [[9]](#References).\nThus with a typical wind of 3-6 m.s$^{-1}$, particles emitted from a source can spread over 15-160 km maximum. Besides, a pixel of satellite data is 3.5 km x 7 km. Knowing this, we make the choice to use a range of 1 km to 100 km characteristic distances and to compute emission factors with different values. "},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"##estimate characteristic distance and surface of a pixel in NO2 pictures\n\nnode = (-67.32431391288841,18.564903861343627)#a corner of the images\nmax_dist_map_long_lat = np.amax(np.linalg.norm(Coords_per_time[4] - node, axis=2, ord = 2))#maximum distance to this corner\nprint(\"The maximum distance in (longitude,latitude) between two points of an image is = {:.4f}\".format(max_dist_map_long_lat))\n\n#characteristic distance\ndist_in_km = 10#characteristic distance we set, in kilometers\nmax_dist_map_in_km = 260#approximation of the maximum distance between two points in kilometers source : google maps\ndist_in_lon_lat = dist_in_km*max_dist_map_long_lat/max_dist_map_in_km#simply use proportions\nprint(\"Characteristic distance within longitude and latitude coordinates change = {:.4f}\".format(dist_in_lon_lat))\n\n#surface of a pixel\nwidth_pixel_long = np.abs(Coords_per_time[0][0,0,0]-Coords_per_time[0][0,1,0])\nwidth_pixel_km = width_pixel_long*max_dist_map_in_km/max_dist_map_long_lat\nheight_pixel_lat = np.abs(Coords_per_time[0][0,0,1]-Coords_per_time[0][1,0,1])\nheight_pixel_km = height_pixel_lat*max_dist_map_in_km/max_dist_map_long_lat\npixel_surface = height_pixel_km*width_pixel_km*pow(10,6)\nprint(\"Width of a pixel (longitude) = {}\".format(width_pixel_long))\nprint(\"Width of a pixel (km) = {:.2f}\".format(width_pixel_km))\nprint(\"Height of a pixel (latitude) = {}\".format(height_pixel_lat))\nprint(\"Height of a pixel (km) = {:.2f}\".format(height_pixel_km))\nprint(\"Surface of a pixel = {:.2f} m^2\".format(pixel_surface))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's select emissions with a gaussian distribution centered on a cluster "},{"metadata":{},"cell_type":"markdown","source":"The gaussian filter is motivated by a physical background: gaussian fields are commonly used to model diffusion processes, which is the case here with a local source and wind creating a flow of molecules."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def Field_evolution_in_a_location(loc_long,loc_lat,field_over_time,coordinates_over_time, plot_map = False):\n '''We use distance in coordinates as a weight to count emissions and store total over time'''\n node = [loc_long,loc_lat]#the center of the filter : source coordinates\n nb_time_steps = field_over_time.shape[0]#amount of time steps\n Weighted_total_over_time = []#we store the sum of tropospheric density pixel after gaussian filtering\n Days_kept = []#to know which days are kept\n for t in range(nb_time_steps):\n dist = np.sum((coordinates_over_time[t] - node)**2, axis=2)#distance^2 between each pixel and the source pixel\n weighted_field = field_over_time[t]*np.exp(-dist/(2*dist_in_lon_lat**2))#we apply weighing\n if(not np.isnan(weighted_field).any()):#check nan values \n Days_kept.append(t)\n if(plot_map):#to plot the 2D fields\n print(t)\n Y,X = np.meshgrid(np.arange(475),np.arange(148))\n plt.pcolormesh(Y,X, weighted_field)\n plt.colorbar()\n plt.show()\n Weighted_total_over_time.append(np.sum(weighted_field, axis = (0,1)))#we append the weighed sum\n return(np.array(Weighted_total_over_time),np.array(Days_kept))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#Let's focus on an example\ntarget_plant_index = 23#example of plant index\ncluster = [23]#example of list of indexes : a cluster of plants\nlong_target, lat_target = float(power_plants[\"longitude\"][target_plant_index]), float(power_plants[\"latitude\"][target_plant_index])\nprint('focus : longitude, latitude = {:.3f}, {:.3f}'.format(long_target,lat_target))\nW_total_per_time, time_kept = Field_evolution_in_a_location(long_target, lat_target,tropospheric_NO2_per_time,Coords_per_time, plot_map = False) ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's create a function to get datetime dates instead of strings"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def to_datetimedates(string_date_list, format_str = '%Y%m%d'):\n '''takes a list of strings and returns a list of datetime dates'''\n Date_list = []\n for date_str_ in string_date_list:\n date_ = datetime.strptime(date_str_, format_str)\n Date_list.append(date_)\n return(np.array(Date_list))\n\ndef to_stringdates(date_list, format_str = '%Y/%m/%d'):\n '''takes a list of datetime dates and returns a list of string dates'''\n Date_string_list = []\n for date_ in date_list:\n date_str_ = date_.strftime(format_str)\n Date_string_list.append(date_str_)\n return(np.array(Date_string_list))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We store the dates of the files without *NaN* in the NO2 dataset:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"dates_NO2 = to_datetimedates(Time_steps,format_str = '%Y/%m/%d' )\ndates_NO2_string = Time_steps\nHours_NO2_instant = to_datetimedates(Hours_NO2_pictures,format_str = '%H%M%S')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def display_date_histogram_per_month(datetime_list):\n '''Takes a list of datetime objects as input and displays a histogram of month repartition of those dates'''\n dtm = lambda x: x.month\n months = list(map(dtm, datetime_list))\n fig, ax = plt.subplots()\n bins = np.arange(1,14)\n ax.hist(months,bins = bins, edgecolor=\"k\", align='left')\n plt.title(\"Repartition of samples considered per month\")\n plt.xlabel(\"Month\")\n plt.ylabel(\"Amount of samples\")\n ax.set_xticks(bins[:-1])\n ax.set_xticklabels([datetime(2020, i, 1).strftime('%b') for i in bins[:-1]])\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"display_date_histogram_per_month(dates_NO2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The following histogram represents the repartition of the time when NO2 pictures were taken (UTC)."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"df = pd.DataFrame({'hour' :Hours_NO2_instant}) \ndf.set_index('hour', drop=False, inplace=True)\ndf = df['hour'].groupby(pd.Grouper(freq='60Min')).count()\nax = df.plot(kind='bar', color='b')\nticklabels = df.index.strftime('%H:%Mh')\nax.xaxis.set_major_formatter(ticker.FixedFormatter(ticklabels))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Remarks on the dates and hours of observed and kept data"},{"metadata":{},"cell_type":"markdown","source":"Most of the acquisition times are between 16 and 17 UTC. This has several consequences:\n\n- The date histogram above highlights the fact that the amount of daily NO2 values we considered after *NaN* processing and interpolation is quite regular (compared to a majority of pictures during summer when excluding any picture with *NaN* values). This enables us to avoid a bias in our results since research papers such as [[12]](#References) have explained that \"[They] expect background influences to be relatively higher at rural than urban sites, and higher in summer (lightning, soil, intercontinental transport; [[13]](#References)) than in winter.\" and that \"There is indeed evidence that free tropospheric NO2 makes a large contribution to OMI NO2 columns and that models underestimate this contribution\".\n

\n- We will need to merge with *GFS* dates and hours in the proper way in next sections (see equation [(4)](#mjx-eqn-convex)).\n

\n- We have to link this instantaneaous observations daily made with the whole NO2 plant emission. We investigate this thorny issue in the following paragraph."},{"metadata":{},"cell_type":"markdown","source":"### C) From instantaneous NO2 to plant emission factor"},{"metadata":{},"cell_type":"markdown","source":"We now need to **derive daily NO2 emissions $e_{plant,day}(n)$ [kg.day$^{-1}$] from a picture of surfacic NO2 taken instantaneously at time $t_{pict}$ for a given plant**. Let $a_{plant}(t)$ [mol] be the amount of NO2 molecules emitted by this plant that are still 'alive' and visible at time $t$ (*i.e.* that are counted in the instantaneaous picture of surfacic density made by the satellite). We can write:\n \n\\begin{equation} \\tag{1}\\label{activity}\n\\frac{[NO]+[NO_2]}{[NO_2]} a_{plant}(t) M_{NO_2} = (e_{plant} * p) (t) = \\int_{\\lambda=0}^{+\\infty} e_{plant}(t-\\lambda)p(\\lambda)\\text{d}\\lambda\n\\end{equation}\n \n- $e_{plant}(t-\\lambda)$ [kg.s$^{-1}$] is the emission rate of NO2 particules emitted a time $t-\\lambda$\n

\n- $p(\\lambda)$ [unitary] is the probability for a particule to live longer than time $\\lambda$.\n

\n- $M_{NO_2}$ [kg.mol$^{-1}$] is the NO2 molar mass.\n

\n- $[NO]\\big/[NO_2] \\simeq 0.32$ is the fraction of NO2 that is chemically converted into NO under urban conditions at noon according to [[14]](#References) and [[3]](#References).\n \nHence the *N02-amount-to-virtual-emitted-mass* conversion factor $k_{mass}$ [kg.mol$^{-1}$] is given by: \n \n$$k_{mass} := M_{NO_2} \\big([N0]+[NO_2]\\big)\\big/[NO_2] \\simeq 6.07 \\times 10^{-2}$$ \n\nLet's write $\\tau$ the expected value lifetime for a particule: \n\n$$\\tau:=\\int_{\\lambda=0}^{+\\infty}\\lambda \\frac{p(\\lambda)-p(\\lambda + \\text{d}\\lambda)}{\\text{d}\\lambda} \\text{d} \\lambda = \\int_{\\lambda=0}^{+\\infty} p(\\lambda)\\text{d} \\lambda$$\n\n**We assume that the probability $p(\\lambda)$ decays quickly comparing to the typical emission rate variation time: $\\frac{1}{e}\\big|\\frac{\\mathrm{d}e}{\\mathrm{d}t}\\big|\\ll \\frac{1}{p}\\big|\\frac{\\mathrm{d}p}{\\mathrm{d}t}\\big|$**. We think that this is quite reasonable since, as we mentioned above, NO2 molecules have a lifetime $\\tau$ that varies from 1 to 8 hours: $\\tau \\sim 2$ hours and one can roughly estimate $p(\\lambda)\\simeq exp(-\\lambda/\\tau)$. This lifetime $\\tau$ is then a lower bound for the resolution of $e_{plant}(t)$ estimation. Hence we can simplify equation [(1)](#mjx-eqn-activity):\n\n$$k_{mass} a_{plant}(t) := \\int_{\\lambda=0}^{+\\infty} e_{plant}(t-\\lambda)p(\\lambda)\\text{d}\\lambda \\simeq e_{plant}(t)\\int_{\\lambda=0}^{+\\infty} p(\\lambda)\\text{d}\\lambda := e_{plant}(t).\\tau $$ \n\nWe also **assume plant emission-rate $e_{plant}(t)$ not to vary over a full day** so that we can extrapolate $e_{plant,day}(n)$ the following way:\n\n$$ e_{plant,day}(n) := \\int_{t \\in day \\, n} e_{plant}(t) \\mathrm{d}t \\simeq T_{day}.e_{plant}(t_{pict}) \\simeq \\frac{T_{day}k_{mass}a_{plant}(t_{pict})}{\\tau} $$\n \n\\begin{equation} \\tag{2}\\label{factor} \\boxed{\ne_{plant,day}(n) \\simeq k_{inst} k_{mass} a_{plant}(t_{pict})}\n\\end{equation}\n \nwhere $k_{inst}:=T_{day}/\\tau \\simeq 12$ is the *instantaneous-to-daily-emission* conversion factor occuring directly in the final result of emission factor in the power-plant contibution, having a heavy dependency $1/\\tau$ on $\\tau$. As $\\tau$ might change according to various conditions (solar exposure, chimical unbalance ...), **the final result for the plant emissions factor $EF_{plant}$ has to be carefully considered with $\\tau$ relative error bar contribution**, see equation [(3)](#mjx-eqn-errorbar) for more details. "},{"metadata":{},"cell_type":"markdown","source":"#### Evolution of the total mass of emissions above the cluster, after *NaN* exclusion"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#constants inherent in NO2 data\nTau = 2 #hours\nprint('Rough NO2 lifetime estimation Tau used : {} hours'.format(Tau))\nk_inst = 24/Tau\nprint('k_inst = {}'.format(k_inst))\nMolar_mass_of_NO2 = 46.0055/1000 #kg/mol\nNO2_to_NOx = 0.32 #chemical balance \nk_mass = (1+NO2_to_NOx)*Molar_mass_of_NO2\nprint('k_mass = {:.2e} kg/mol'.format(k_mass))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#process one particular cluster\nmean_on_kept_period = k_inst*Molar_mass_of_NO2*pixel_surface*np.mean(W_total_per_time)##mean of daily NO2 mass close to the plant\nplt.plot(k_inst*Molar_mass_of_NO2*pixel_surface*W_total_per_time)\nplt.plot(mean_on_kept_period*np.ones(np.shape(W_total_per_time)))\nplt.xlabel(\"Days kept after NaN values processing\")\nplt.ylabel(\"Mass of tropospheric NO2: kg\")\nplt.title(\"Mass of tropospheric vertical column of NO2 \\n linked to daily emissions above the cluster\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Mix with plant dataset to get an approximation of the emission factor of a cluster"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"print(\"Mean of mass of tropospheric NO2 linked to daily emissions above the cluster: {:.2f} kg\".format(mean_on_kept_period))\nfirst_approx_emission = 365*(1+NO2_to_NOx)*mean_on_kept_period\ngeneration_per_plant_cluster = power_plants[\"estimated_generation_gwh\"][cluster]\ngeneration_cluster = np.sum(np.array(generation_per_plant_cluster))\nfirst_approx_emission_factor = first_approx_emission/generation_cluster\nprint(\"NOx emission (average) over one year linked to the cluster: {:.2f} kg\".format(first_approx_emission))\nprint(\"NOx emission factor (average) over one year linked to the cluster: {:.2f} NOx kg/GWh\".format(first_approx_emission_factor))\nprint(\"Documented values are between 100 (300 in the introduction table) \\n and 1000 (4000 in the same table) kg/GWh\") \n#see other documented values here : https://books.google.fr/books?id=2L-OIrTBiRMC&pg=PA176&lpg=PA176&dq=NO2x+kg/GWh&source=bl&ots=HpOZw0yPip&sig=ACfU3U3i0w4Ya5ex2vfwOwjNMk27JiG7qg&hl=fr&sa=X&ved=2ahUKEwjws4H_vOPnAhVoDWMBHds2BdAQ6AEwAHoECAcQAQ#v=onepage&q=NO2x%20kg%2FGWh&f=false","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### D) Extend approach to clusters and sum to get a first estimation for the whole island"},{"metadata":{},"cell_type":"markdown","source":"This section is a generalization of the previous one: we create a set of clusters grouping fossil plants regarding their coordinates and the characteristic distance used. We compute emissions for each cluster and get elecricity generation too. In the end, we combine both information to estimate emission factors."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#created by hand with the map\nClusters = [[16],[2,12,14],[23],[18],[17,3,30],[13],[11],[19],[28,15,10,1]]#one cluster : one list of indexes","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def barycenter(coord_list):\n '''gets a list of coordinates and returns the coordinates of the barycenter'''\n return(np.sum(np.array(coord_list), axis = 0)/len(coord_list))\n\ndef index_to_coordinates(ind):\n '''gives the coordinates in longitude, latitude of the plant with index = ind in the dataframe'''\n long_,lat_ = float(power_plants[\"longitude\"][ind]), float(power_plants[\"latitude\"][ind])\n return([long_,lat_])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's store the barycentre coordinates of each cluster:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Barycenters = []\nfor clust in Clusters:\n coord_clust = [index_to_coordinates(ind) for ind in clust]\n bar = barycenter(coord_clust)\n Barycenters.append(bar)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We define a function to generalize the approach described above: for each cluster we compute the evolution of NO2 density per day, weighed by the gaussian filter and computed as a total (sum) in the troposphere."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def track_NO2_emission_density_per_cluster_per_time(Clusters_,Barycenters_):\n '''Takes a list of cluster indexes and the coordinates of the barycenter of each cluster, returns the total NO2 emissions'''\n NO2_per_cluster_per_time = []\n for k in range(len(Clusters_)):\n print(\"Dealing with cluster : {}\".format(k))\n barycenter_ = Barycenters_[k]\n W_total_per_time, _ = Field_evolution_in_a_location(barycenter_[0],barycenter_[1],tropospheric_NO2_per_time,Coords_per_time)\n NO2_per_cluster_per_time.append(W_total_per_time)\n return(np.array(NO2_per_cluster_per_time))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's apply it to our clusters:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Time_series_NO2_density_per_cluster = track_NO2_emission_density_per_cluster_per_time(Clusters, Barycenters)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's see evolution over time and distribution of emissions"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#evolution over time per cluster\nMean_per_cluster = []\nfor clust_ind in range(len(Clusters)):\n mean_on_kept_period = k_inst*Molar_mass_of_NO2*pixel_surface*np.mean(Time_series_NO2_density_per_cluster[clust_ind])##mean of daily NO2 mass close to the plant\n Mean_per_cluster.append(mean_on_kept_period)\n plt.plot(k_inst*Molar_mass_of_NO2*pixel_surface*Time_series_NO2_density_per_cluster[clust_ind])\n plt.plot(mean_on_kept_period*np.ones(np.shape(Time_series_NO2_density_per_cluster[clust_ind])))\n plt.xlabel(\"Days kept after NaN values processing\")\n plt.ylabel(\"Mass of tropospheric NO2: kg\")\n plt.title(\"Mass of tropospheric vertical column of NO2 \\n linked to daily emissions, cluster : {}\".format(clust_ind))\n plt.show()\n \nMean_per_cluster = np.array(Mean_per_cluster)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#distribution per cluster\nfor clust_ind in range(len(Clusters)):\n sns.distplot(k_inst*Molar_mass_of_NO2*pixel_surface*Time_series_NO2_density_per_cluster[clust_ind])\n plt.xlabel(\"Mass of troposheric NO2 (kg)\")\n plt.ylabel(\"Density function\")\n\nplt.title(\"Distribution of mass of tropospheric vertical column of NO2 \\n linked to daily emissions per cluster\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Mix with plant dataset to get an approximation of the emission factor"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#emissions\nfirst_approx_emission_per_clust = 365*(1+NO2_to_NOx)*Mean_per_cluster\nsns.distplot(first_approx_emission_per_clust/1000)\nplt.title(\"Distribution of NOx emission per cluster\")\nplt.xlabel(\"Mass of NOx : tons\")\nplt.ylabel(\"Density\")\nplt.show()\n\n#electricity generation\ngeneration_per_clust = []\nfor cluster in Clusters:\n generation_inside_clust = np.array(power_plants[\"estimated_generation_gwh\"][cluster])\n generation_per_clust.append(np.sum(generation_inside_clust)) \ngeneration_per_clust = np.array(generation_per_clust)\nprint(\"Annual total production per cluster \\n {} GWh\".format(np.array(generation_per_clust,dtype=np.float32)))\n\nsns.distplot(generation_per_clust)\nplt.title(\"Distribution of electricity generation per cluster\")\nplt.xlabel(\"Electricity generated : GWh\")\nplt.ylabel(\"Density\")\nplt.show()\n\n#emission factor\nfirst_approx_emission_factor_global = first_approx_emission_per_clust/generation_per_clust\n#print(\"Mean NOx emission factor over one year per cluster : = {} NOx kg/GWh\".format(first_approx_emission_factor_global))\nprint(\"Mean per cluster = {:.2f} NOx kg/GWh\".format(np.mean(first_approx_emission_factor_global)))\nprint(\"Documented values are between 100 (300 in the introduction table) and 1000 (4000 in the same table) kg/GWh\")\n#see for other values : https://books.google.fr/books?id=2L-OIrTBiRMC&pg=PA176&lpg=PA176&dq=NO2x+kg/GWh&source=bl&ots=HpOZw0yPip&sig=ACfU3U3i0w4Ya5ex2vfwOwjNMk27JiG7qg&hl=fr&sa=X&ved=2ahUKEwjws4H_vOPnAhVoDWMBHds2BdAQ6AEwAHoECAcQAQ#v=onepage&q=NO2x%20kg%2FGWh&f=false\n\nsns.distplot(first_approx_emission_factor_global)\nplt.title(\"Distribution of emission factor per cluster\")\nplt.xlabel(\"Emission factor : NOx kg/GWh\")\nplt.ylabel(\"Density\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's deduce an emission factor for the whole island:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"sum_generation = power_plants[\"estimated_generation_gwh\"].sum()\nprint(\"Mean emission factor for electricity generation in Puerto Rico = {:.0f} NOx kg/GWh\".format(first_approx_emission_per_clust.sum()/sum_generation))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Comments on those results and error bar estimation"},{"metadata":{},"cell_type":"markdown","source":"As explained in the previous remarks, this first rough estimation for emission factor was carried out using a lot of hypotheses and approximations whose effects have to bo taken into account in this result interpretation. \n\n- First point : on average and in total, **we tend to overestimate emissions and thus emissions factors with this gaussian approach**. This can be linked to various assumptions and factors:\n\n - As we mentioned previously, we must have in mind that natural background could account for a large part of tropospheric NO2 detected. We will investigate this point in the next section.\n \n - One can notice with a rapid view over the map that in most cases, plants are surrounded by cities, or even by aiports. It seems difficult to distinguish those sources of NOx emissions with the current model. A way to address this issue could be to focus only on isolated plants and to select a small characteristic distance. Nonetheless, a compromise has to be taken not to underestimate emissions because of a too small distance, and generalization with only a small sample of the dataset can be questioned.\n \n - To extrapolate full daily emissions from instantaneous pictures of NO2, we multiplied instantaneous NO2 density with a *instantaneous-to-daily-emission* factor $k_{inst}:=T_{day}/{\\tau}$ (see equation [(2)](#mjx-eqn-factor)), based on the questionable hypothesis of a constant emission for each power-plant over a full day, and accurate knowledge of average NO2 lifetime $\\tau$. Nonetheless one could imagine that power-plants are working less during the night. This is thus an upper bound to plant emissions.\n \n \n- We can estimate the relative errobar of the emission factor $EF_{cl,y}$ of each cluster $cl$ averaged over a year using days set $Days$ of size $N_{days}$ defined as:\n\n \\begin{equation*}\n EF_{cl,y} = \\frac{e_{cl,y}}{prod_{cl,y}} = \\frac{365}{N_{days}}\\frac{\\sum_{n \\in Days} e_{cl,day}(n)}{prod_{cl,y}}\n \\end{equation*}\n \n Since the production data $prod_{cl,y}$ is supposed to be reliable, we only take into account the pollution relative uncertainty that basically depends on the 3 factors $k_{inst}$, $k_{mass}$ and $a(t)$ (as seen in equation [(2)](#mjx-eqn-factor)). Because we set once and for all the two first factors $k_{inst}$ and $k_{mass}$, we introduce here a constant bias that cannot be reduced using a statistical averaging. The same issue comes for $a(t)$ that tends to be constantly biased because of positive background contribution:\n \n \\begin{equation} \\label{errorbar} \\tag{3}\n \\frac{\\Delta EF_{cl,y}}{EF_{cl,y}} = \\frac{\\Delta e_{cl,y}}{e_{cl,y}} = \\sqrt{\\bigg(\\frac{\\Delta \\tau}{\\tau}\\bigg)^2+ \\bigg(\\frac{\\Delta k_{mass}}{k_{mass}} \\bigg)^2 + \\bigg(\\frac{\\Delta a}{a}\\bigg)^2}\n \\end{equation}\n \n where $\\big(\\Delta \\tau / \\tau \\big)^2$, $\\big(\\Delta k_{mass} / k_{mass} \\big)^2$ and $\\big(\\Delta a / a \\big)^2$ are averaged over the $Days$ sample. We can make some crude estimations in order to assess our result errobar:\n \n \\begin{array}{|c|c|c|}\n \\hline\n \\Delta \\tau / \\tau \\sim 100\\% & \\Delta k_{mass} / k_{mass} \\sim 20\\% & \\Delta a / a \\sim 150\\% \\\\\n \\hline\n \\end{array}\n $$ $$\n $$ \\boxed{\\frac{\\Delta EF_{cl,y}}{EF_{cl,y}} \\sim 180 \\% } $$\n \n \n- We have had to exclude or to modify some values : the estimated generation of the coal and hydro plants were not consistent and we had a lot of *NaN* values in satellite data (maybe due to cloud coverage that entails some reflection, since NO2 density is deemed thanks to optic properties).\n\n\n- Clusters have to be thought in accordance with the characteristic distance. Indeed, if this latter is greater or almost equal to the distance between two clusters, we will face some issues counting twice some emissions. Clusters must thus be modified having in mind this principle when the characteristic distance is altered.\n\n\n- According to dataset description, the estimation of generated electricity in the power plants file is about 2015, and our NO2 satellite data is linked to 2018. As an approximation, we assumed that emissions data for 2015 and 2018 were close. Some ways to deal with this date difference would be to incorporate a new dataset for 2018, or to use activity models over time to deduce from 2015 data an approximation of 2018 data."},{"metadata":{},"cell_type":"markdown","source":"## IV - Second attempt: prepare a broader dataset with horizontal wind, temperature and relative humidity"},{"metadata":{},"cell_type":"markdown","source":"*Go back to table of contents [⤴](#Content-of-the-rest-of-this-notebook)*"},{"metadata":{},"cell_type":"markdown","source":"The aim of this section is to deal with limitations of the previous one we noticed above. \n\nIndeed, as depicted in [[15]](#References), humidity and temperature are features that play no mean part in natural NO2 emissions. The importance of this background emission was mentioned previsouly. \n\n\nBesides, our model with gaussian masks was underpinned by the assumption of a typical wind speed of 3-6 m.s$^{-1}$, and did not take into account this variable to compute emissions from NO2 2D field. \n\nWe are thus trying here to improve the previous work, taking into account wind speed and orientation, relative humidity and temperature as additional features. We fit 1D models taking into account those new variables, as bolstered by [[9]](#References). This latter was inspiring as a valuable source of information."},{"metadata":{},"cell_type":"markdown","source":"### A) Extract *GFS* 2D fields and match the dates with NO2"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/'\n\ngfs_coords_per_time = {12 : [], 18 : []}\ngfs_wind_per_time = {12 : [], 18 : []}\ngfs_temperature_per_time = {12 : [], 18 : []}\ngfs_humidity_per_time = {12 : [], 18 : []}\ngfs_time_laps = {12 : [], 18 : []}\n\nfor extract_iteration, filename in enumerate(os.listdir(path)):\n date_ = filename[4:14]\n if date_[:-2] in Time_laps[:,0]: #we extract only on the days of interest\n hour_ = int(date_[-2:])\n if hour_ in [12, 18]: #we extract only hours close to NO2 acquisition hour \n #print('Extracting file number {} : '.format(extract_iteration)+filename)\n gfs_time_laps[hour_].append(date_)\n path_ = path+filename\n gfs_coords_per_time[hour_].append(coordinates(path_))\n gfs_wind_per_time[hour_].append(rio.open(path_).read()[3:5])\n gfs_temperature_per_time[hour_].append(rio.open(path_).read()[0])\n gfs_humidity_per_time[hour_].append(rio.open(path_).read()[2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sort data by date:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"gfs_dates = dict()\ngfs_timestamp = dict()\n\nfor hour in [12, 18]: \n gfs_dates[hour] = to_datetimedates(np.array(gfs_time_laps[hour]),\n format_str='%Y%m%d%H')\n gfs_timestamp[hour] = np.array([datetime.timestamp(date)\n for date in gfs_dates[hour]])\n gfs_coords_per_time[hour] = np.array(gfs_coords_per_time[hour])\n gfs_wind_per_time[hour] = np.array(gfs_wind_per_time[hour])\n gfs_temperature_per_time[hour] = np.array(gfs_temperature_per_time[hour])\n gfs_humidity_per_time[hour] = np.array(gfs_humidity_per_time[hour])\n\n print('sorting the data for hour '+str(hour))\n total_data = {instant: (date, coords, wind, temperature, humidity)\n for instant, date, coords, wind, temperature, humidity\n in zip(gfs_timestamp[hour],\n gfs_dates[hour],\n gfs_coords_per_time[hour],\n gfs_wind_per_time[hour],\n gfs_temperature_per_time[hour],\n gfs_humidity_per_time[hour])}\n \n gfs_timestamp[hour] = np.array(sorted(gfs_timestamp[hour]))\n gfs_dates[hour] = np.array([total_data[instant][0]\n for instant in gfs_timestamp[hour]])\n gfs_coords_per_time[hour] = np.array([total_data[instant][1]\n for instant in gfs_timestamp[hour]])\n gfs_wind_per_time[hour] = np.array([total_data[instant][2]\n for instant in gfs_timestamp[hour]])\n gfs_temperature_per_time[hour] = np.array([total_data[instant][3]\n for instant in gfs_timestamp[hour]])\n gfs_humidity_per_time[hour] = np.array([total_data[instant][4]\n for instant in gfs_timestamp[hour]])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check date information about those files"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"if (gfs_dates[12].shape == np.unique(gfs_dates[12]).shape) and \\\n (gfs_dates[18].shape == np.unique(gfs_dates[18]).shape):\n print(\"No duplicates in GFS data\")\nelse:\n print(\"Duplicates, need to process dates\")\n \nif gfs_dates[12].shape == gfs_dates[18].shape:\n print(\"Same amount of data for 12 UTC and 18 UTC\")\n if len(gfs_dates[12]) == len(timestamp):\n print(\"Data is known for day of NO2 data\")\n else:\n print(\"Data misses for some day of NO2 data\")\nelse:\n print(\"Issue : different amount of data for 12 UTC and 18 UTC\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's process the wind vector field stored\n\nWe know that the U component is related to East-West axis and oriented positive towards East, V is related to North-South axis and oriented positive towards North. Besides, satellite pictures which are provided are oriented East-West and North-South too."},{"metadata":{},"cell_type":"markdown","source":"We need first to compute the norm of the wind field, to know whether the previous hypothesis on diffusion was consistent."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def plot_wind_norm(wind_field):\n '''Takes a wind 2D field as input and displays its scalar norm field through a color map'''\n norm_field = np.linalg.norm(wind_field, axis = 0, ord = 2)\n Y,X = np.meshgrid(np.arange(475),np.arange(148))\n plt.pcolormesh(Y,X, norm_field)\n plt.xlabel('Pixel index in East-West direction')\n plt.ylabel('Pixel index in North-South direction')\n plt.title('Wind norm over the island')\n cb = plt.colorbar()\n cb.set_label('m/s')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can display the scalar field over the island for the first wind file:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"plot_wind_norm(gfs_wind_per_time[12][0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that for some particular days and locations, the speed is higher than we assumed. Let's compute an evolution over time of the average norm over the island:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def spatial_average_over_time(field, vector_field = True):\n '''Takes a field over time above Puerto Rico as input, and returns its spatial average over time'''\n '''If the field is a vector field, then its average norm is returned over time (vector_field = True)'''\n if(vector_field):\n return(np.mean(np.linalg.norm(field, axis = 1, ord = 2), axis = (1,2)))\n else:\n return(np.mean(field, axis = (2,3)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"spatial_wind_norm_average_over_time = {12 : spatial_average_over_time(gfs_wind_per_time[12]),\n 18 : spatial_average_over_time(gfs_wind_per_time[18])}","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"plt.plot(spatial_wind_norm_average_over_time[18])\nplt.xlabel(\"Time index in the GFS dataset\")\nplt.ylabel(\"Spatial average of wind norm (m/s)\")\nplt.title(\"Evolution of the spatial average of the wind norm above the island at 18 UTC\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"sns.distplot(spatial_wind_norm_average_over_time[18])\nplt.title(\"Distribution of spatial average of the wind norm at 18 UTC above the island over one year\")\nplt.xlabel(\"Spatial average of wind norm (m/s)\")\nplt.ylabel(\"Density\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"print(\"Spacetime average of the wind norm above the island at 12 UTC = {:.4f} m/s\".format(np.mean(spatial_wind_norm_average_over_time[12])))\nprint(\"Spacetime average of the wind norm above the island at 18 UTC = {:.4f} m/s\".format(np.mean(spatial_wind_norm_average_over_time[18])))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Though the distribution is not gaussian, we can notice that the evolution of the spatial average of the wind norm and its distribution enable to confirm the assumption taken in the previous section to select the characteristic distance and to compute the emission factor. Our first approximation was not absurd."},{"metadata":{},"cell_type":"markdown","source":"#### We are going to deal with streamlines before creating a model"},{"metadata":{},"cell_type":"markdown","source":"First, let's deal with wind orientation and display stream lines, that are tangent to the local wind."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def plot_wind_streamline(wind_field):\n '''Takes a 2D wind field as input and returns its streamlines'''\n Y,X = np.meshgrid(np.arange(475),np.arange(147,-1,-1))\n U, V = wind_field[0], wind_field[1]\n plt.title(\"Streamlines of the wind above the island\")\n plt.xlabel('Pixel index in East-West direction')\n plt.ylabel('Pixel index in North-South direction')\n strm = plt.streamplot(Y, X, U, V, color=np.sqrt(U**2+V**2),linewidth= 2, cmap=plt.cm.viridis)\n cb = plt.colorbar(strm.lines)\n cb.set_label('m/s')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are going to display the 2 streamline pictures of the wind inherent in one day of wind data. \n\nOur concern is the following : if wind direction evolves too much over a day, we won't be easily able to use wind direction to process NO2 daily data taken between 15h and 18h UTC."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"day_index = 1 #day index in wind observation\nfor hour in [12, 18]:\n print(\"We display streamlines for date and time : {}\".format(gfs_dates[hour][day_index]))\n plot_wind_streamline(gfs_wind_per_time[hour][day_index])\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are going to merge the two data sets of *GFS* given for 12 UTC and 18 UTC for each day in order to have an approximation of the wind, temperature and humidity during the intantaneous NO2 acquistion made at time $t_{pict}$ stored in *Hours_NO2_instant*. Every field $f$ will be interpolated as follow: \n\n\\begin{equation}\\label{convex}\\tag{4} \nf(x,y,t_{pict}) \\simeq \\big(1-\\alpha(t_{pict})\\big) f(x,y,t_0) + \\alpha(t_{pict}) f(x,y,t_1)\n\\end{equation}\n\n$$\\text{where } \\alpha(t_{pict})=\\frac{t_{pict}-t_0}{t_1-t_0} $$\n\nfor $t_{pict}$ (given in hours UTC) in the range : $[t_0=12h, t_1=18h]$"},{"metadata":{},"cell_type":"markdown","source":"We can see first that the *GFS* coordinates do not depend on time, so they can be kept in a single array *gfs_coords*:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"gfs_coords = gfs_coords_per_time[12][0]\nprint(np.max(abs(gfs_coords_per_time[12]-gfs_coords)))\nprint(np.max(abs(gfs_coords_per_time[18]-gfs_coords)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So we can proceed to direct convex combinations over wind, temperature and humidity arrays:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"alphas = np.array(list(map(lambda hour :\n (hour.hour+hour.minute/60+hour.second/3600-12)/(18-12),\n Hours_NO2_instant)))\nif np.min(alphas)<0 or np.max(alphas)>1:\n print('Issue : there is an acquisition time not between 12h and 18h UTC')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"wind_per_day = np.array(list(map(lambda x : x[0]*x[1],\n zip(1-alphas, gfs_wind_per_time[12]))))+\\\n np.array(list(map(lambda x : x[0]*x[1],\n zip(alphas, gfs_wind_per_time[18]))))\ntemperature_per_day = np.array(list(map(lambda x : x[0]*x[1],\n zip(1-alphas, gfs_temperature_per_time[12]))))+\\\n np.array(list(map(lambda x : x[0]*x[1],\n zip(alphas, gfs_temperature_per_time[18]))))\nhumidity_per_day = np.array(list(map(lambda x : x[0]*x[1],\n zip(1-alphas, gfs_humidity_per_time[12]))))+\\\n np.array(list(map(lambda x : x[0]*x[1],\n zip(alphas, gfs_humidity_per_time[18]))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"day_index = 0\nprint(\"Estimated streamlines during NO2 acquisition for date {} and time {} UTC\".format(\\\n Time_steps[day_index],Hours_NO2_instant[day_index].strftime('%Hh:%M')))\nplot_wind_streamline(wind_per_day[day_index])\nplt.show()\nplt.title(\"Corresponding temperature (°C)\")\nplt.imshow(temperature_per_day[day_index])\nplt.colorbar()\nplt.show()\nplt.title(\"Corresponding relative humidity (%)\")\nplt.imshow(humidity_per_day[day_index])\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### B) From 2D fields to 8 1D fields"},{"metadata":{},"cell_type":"markdown","source":"#### We create functions to extract from a 2D field some 1D fields along specific directions and in a perimeter of a location"},{"metadata":{},"cell_type":"markdown","source":"### Main ideas and motivations"},{"metadata":{},"cell_type":"markdown","source":"Let's define the aim and the motivations of the next functions:\n\nOur aim is to compress our very granular 2D fields into few smoother one dimensional arrays to fit physical models. Several aspects are bound to compete and we have to take them into account before we chose the way we extract information:\n\n- As you could see on the previous plotted streamlines of the wind above the island, they are not uniform in space neither time stationary: between 12 UTC and 18 UTC you can see some variations. Therefore, our convex combination [4](#mjx-eqn-convex) to get the wind field at the NO2 acquisition time is limited in terms of accuracy. We believe we can get an idea of the discrete direction in a range of cardinal points, but the **continous value has a limited sense**.\n\n\n- In addition, we need to think about the way we extract the NO2 once we have an estimation of the local wind direction at the time of NO2 acquisition.\n - One can think about a rectangle along the axis of this latter direction. The issue is that if we do not estimate well that direction, or if wind is not really stationary, we would select the wrong information: we would consider pixels that do not correspond to power-plant-related emissions (moving away from the barycentre of the cluster some pixels would be out of the rectangle, and very close to the barycentre we would select too many pixels including some having no link with power-plants). \n - We thus thought about another approach : to select the NO2 pixels in **an eighth of a disc** along the discrete cardinal point that is the closest to wind direction. In this way, we accept a certain lack of precision on the direction and we preferably select NO2 based on barycentre coordinates and wind information. This choice will entail some difficulties when it comes to compute integrals, but seems to respect more the physical description of the diffusion problem. Indeed, in case of weak wind or 'isotropic' wind (*i.e.* without any main direction standing out during the diffusion), the problem appears to be symmetric with respect to the plant localisation and **encourages us to use a cylindrical description centered over the plant** (more broadly the cluster barycentre). Furthermore, when having a main wind direction, the diffusion is mainly gathered in a 2D cone starting from the plant, not a rectangle.\n\n\n- Another point to be questioned now is the radial size and resolution we set for this eighth of a disc. The radius will be set to the maximum distance of a pixel from the barycentre of the cluster we believe to be power-plant related. To decide on this topic, we have in mind the average wind speed 3-4 m.s$^{-1}$ and the average lifetime $\\tau \\sim 2$ hours of NO2 molecules that entail a range of 10-30 km. **We have set a radius of 15 km**. Because of the extreme granularity of 2D fields (3.5km x 7km for the NO2) over such a small cone, **we discretized it only in 10 evenly spaced radiuses**.\n\n\n- Our choice of the cylindric approach entails a compromise. Indeed, though it seems more natural and a better way to get the proper information, integration is not straightforward. We will compute 1D lineic densities from 2D surfacic densities, by averaging in the orthogonal direction to the principal wind direction. To do so, we thought we had two options :\n - Compute averages and integrals rigorously on every portion of the eighth of a disc, taking into account that since we separate with regular intervals along this principal direction, the portions we consider increase when we move apart from the barycentre of the cluster. This process induces a kind of distortion that would have to be taken into account in physical models below: the surfacic density averaged and multiplied by the orthogonal distance of the portion considered to get the linear density would be increased away from the barycentre and would be tightened close to this barycentre. We thought correction in the physical model would not be easily implemented. We thus did not chose this approach.\n - Compute surfacic densities over those unequal portions of eighth of a disc, averaging on increasing surfaces. Nonetheless, when it comes to compute a lineic density, we correct the distortion induced by the increase of the orthogonal distance with respect to the principal direction when we move away from the barycentre by multiplying by a dimensionless factor. This strategy can be thought geometrically: though we select information over an eighth of a disc, we compute lineic density as if we it were a rectangle. Thus we integrate in the orthogonal direction to the principal wind direction by mutliplying those average surfacic values by a constant distance (we select this latter as the mean orthogonal distance of the eighth of a disc we consider). We believe this approximation is a way to select plant-related emissions in a compromise to get data that can be fitted by traditional physics models."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Convention_directions = np.array(['North', 'South', 'East', 'West', 'N-E', 'N-W', 'S-E', 'S-W'])\nConvention_vectors = np.array([(0,1),(0,-1),(1,0),(-1,0),(1/np.sqrt(2),1/np.sqrt(2)),(-1/np.sqrt(2),1/np.sqrt(2)),(1/np.sqrt(2),-1/np.sqrt(2)),(-1/np.sqrt(2),-1/np.sqrt(2))])\n\ndef get_surrounding(center_coords_, radius_, field_values_, field_coordinates_, plot_mask_ = False):\n '''Extracts the surrounding of a point with a radius, and get field values in this surrounding'''\n '''Radius : km, field_coordinates and center_cords : (lon,lat), field_values : 2D field'''\n Radius_lon_lat = radius_*max_dist_map_long_lat/max_dist_map_in_km#conversion as mentioned in the previous section\n Disk_mask = (np.linalg.norm(field_coordinates_ - center_coords_, axis=2, ord = 2) <= Radius_lon_lat).astype(int) \n if(plot_mask_):\n if(field_coordinates_[0,0][1]>field_coordinates_[1,0][1]):\n Y,X = np.meshgrid(np.arange(475),np.arange(147,-1,-1))\n else:\n Y,X = np.meshgrid(np.arange(475),np.arange(148))\n plt.pcolormesh(Y,X, Disk_mask*field_values_)\n plt.axis('equals')\n plt.colorbar()\n return Disk_mask*field_values_\n\ndef get_8_directions(center_coords,radius,field_values,field_coordinates, \n plot_extract=False, nb_intervals=10):\n '''Computes 8 density lines of the 2D scalar field around the point with coordinates : center_coords and radius limits'''\n '''Lines : towards North, South, East, West, N-E, N-W, S-E, S-W using (longitude, latitude) coordinates'''\n \n translated_coordinates = field_coordinates-center_coords#origin : center_coords\n surronding_values = get_surrounding(center_coords,radius,field_values,field_coordinates)\n N_vector = (0,1)\n S_vector = (0,-1)\n E_vector = (1,0)\n W_vector = (-1,0)\n NE_vector = (1/np.sqrt(2),1/np.sqrt(2))\n NW_vector = (-1/np.sqrt(2),1/np.sqrt(2))\n SE_vector = (1/np.sqrt(2),-1/np.sqrt(2))\n SW_vector = (-1/np.sqrt(2),-1/np.sqrt(2))\n directions = [N_vector,S_vector,E_vector,W_vector,NE_vector,NW_vector,SE_vector,SW_vector]\n scalar_products = np.array([np.dot(translated_coordinates,direction) for direction in directions])\n Directional_densities = []\n \n for direction_index in range(len(directions)):\n \n #we select scalar products and values in the direction and below radius\n mask = (np.argmax(scalar_products, axis = 0) == direction_index).copy()\n selected_locations = (scalar_products[direction_index][mask])[surronding_values[mask] > 0]\n selected_values = (surronding_values[mask])[surronding_values[mask] > 0]\n \n #we create a discrete line to compute a density\n discrete_intervals = np.linspace(0,radius*max_dist_map_long_lat/max_dist_map_in_km, nb_intervals)\n interval_correspondance = np.digitize(selected_locations,discrete_intervals)\n line_data = pd.DataFrame(columns = [\"Values\", \"interval_id\"], data = np.transpose(np.array([selected_values,interval_correspondance])))\n interval_data = pd.DataFrame(columns = [\"interval_index\", \"upper_bound\"], data = np.transpose(np.array([np.arange(nb_intervals),discrete_intervals])))\n discretize_density = line_data.groupby([\"interval_id\"]).mean()#we select the mean per interval along the direction\n density_with_metric = discretize_density.join(interval_data.set_index(\"interval_index\"))\n Directional_densities.append(density_with_metric)\n \n if(plot_extract):\n Y,X = np.meshgrid(np.arange(475),np.arange(148))\n plt.title(\"Conical extraction in {} direction\".format(Convention_directions[direction_index]))\n plt.pcolormesh(Y,X, mask.astype(int)*surronding_values)\n plt.show()\n plt.plot(density_with_metric[\"upper_bound\"],density_with_metric[\"Values\"])\n plt.title(\"1D density result in {} direction\".format(Convention_directions[direction_index]))\n plt.show()\n return Directional_densities\n\ndef vect_to_direction(vect):\n '''Gives the index of the direction corresponding to the vector in North, South, East, West, N-E, N-W, S-E, S-W conventions'''\n N_vector = (0,1)\n S_vector = (0,-1)\n E_vector = (1,0)\n W_vector = (-1,0)\n NE_vector = (1/np.sqrt(2),1/np.sqrt(2))\n NW_vector = (-1/np.sqrt(2),1/np.sqrt(2))\n SE_vector = (1/np.sqrt(2),-1/np.sqrt(2))\n SW_vector = (-1/np.sqrt(2),-1/np.sqrt(2))\n directions = [N_vector,S_vector,E_vector,W_vector,NE_vector,NW_vector,SE_vector,SW_vector]\n scalar_products = np.array([np.dot(vect,direction) for direction in directions])\n return(np.argmax(scalar_products))\n\ndef get_opposite_direction(ind_dir):\n \"\"\"takes as input an index corresponding to one direction in the Convention_directions array\"\"\"\n \"\"\"Returns the index of the opposite direction in the same array\"\"\"\n Opposites = [1,0,3,2,7,6,5,4]#South is index 1 and opposite to North which is index 0 for instance\n return Opposites[ind_dir]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Now we are going to study linkages between one-directional densities and wind direction"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"day_index = 39#example of day index in the list of dates\ncluster_index = 4#cluster index\nrad = 15#perimeter : radius of the disk we consider\n\nprint(\"We focus on day : {}\".format(Time_steps[day_index]))\n\nNO2_test = get_surrounding(Barycenters[cluster_index],\n rad,\n tropospheric_NO2_per_time[day_index],\n Coords_per_time[day_index])\ndensities_NO2_1D = get_8_directions(Barycenters[cluster_index],\n rad,\n tropospheric_NO2_per_time[day_index],\n Coords_per_time[day_index],\n plot_extract = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"for dir_index in range(len(densities_NO2_1D)):\n density_dir = densities_NO2_1D[dir_index]\n plt.plot(density_dir[\"upper_bound\"]*max_dist_map_in_km/max_dist_map_long_lat,\n density_dir[\"Values\"]*Molar_mass_of_NO2*pixel_surface)\nplt.legend(Convention_directions,bbox_to_anchor=(1, 1))\nplt.xlabel(\"Distance from the cluster barycenter: km\")\nplt.ylabel(\"Mass of tropospheric NO2: kg\")\nplt.title('NO2 density above a cluster {} for date {} and time {} UTC'.format(cluster_index,\n Time_steps[day_index],\n Hours_NO2_instant[day_index].strftime('%Hh:%M')))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"wind_distribution = {direction : 0 for direction in Convention_directions}\n\nfor day_index in range(n_t):\n wind_on_location = get_surrounding(Barycenters[cluster_index],\n rad,\n wind_per_day[day_index],\n gfs_coords)\n wind_on_location = wind_on_location.reshape(2,-1)\n wind_on_location = wind_on_location[:,np.linalg.norm(wind_on_location,\n axis = 0,\n ord = 2)>pow(10,-2)]\n mean_wind_on_location = np.mean(wind_on_location, axis = 1)\n dir_index = vect_to_direction(mean_wind_on_location)\n wind_distribution[Convention_directions[dir_index]] += 1\n #print(\"Wind direction for date {} and time {} UTC : {}\".format(Time_steps[day_index],\n # Hours_NO2_instant[day_index].strftime('%Hh:%M'),\n # Convention_directions[dir_index]))\n\nplt.bar(wind_distribution.keys(), wind_distribution.values())\nplt.title(\"Distribution over the year of the main wind direction\", fontsize=18, pad=15)\nplt.xlabel(\"Wind main direction during NO2 acquisition\", fontsize=15)\nplt.ylabel(\"Number of days\", fontsize=15)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the main wind direction above power-plants is West, with a huge majority. In the plot above where NO2 density is displayed per discrete direction around a specific cluster, one can notice that the West 1D density is higher and has a specific form: we can see a kind of diffusion.\nThis impression is formalised just below."},{"metadata":{},"cell_type":"markdown","source":"### Description of the next approach"},{"metadata":{},"cell_type":"markdown","source":"We can notice that the density of tropospheric NO2 evolves in a different way along the wind direction compared to other directions. With functions defined above we have now the opportunity to create a dataframe we will complete with temperature and humidity information.\n\nFor each cluster of power-plants and each day in cleaned NO2 data, we process in the following way:\n$$ $$\n - We use cylindrics coordinates $(r,\\theta)$ with respect to its barycentre.\n$$ $$\n - We compute the average wind $\\vec{v}(t)$ in a specific surrounding of the cluster of power-plants.\n$$ $$ \n - We find the closest discrete direction $\\vec{u} \\in \\big\\{N,NW,W,SW,S,SE,E,NE \\big\\}$ that maximizes $\\vec{v}(t) \\cdot \\vec{u}$. Thus we select the closest discrete direction to the wind direction.\n$$ $$ \n - We store the mere scalar $\\vec{v}(t) \\cdot \\vec{u}$ to account for wind contribtion in the model.\n$$ $$ \n - We compute the **lineic density** $NO2_{\\vec{u}}(r,t)$ in the chosen direction $\\vec{u}$, temperature above ground and relative humidity above ground **over an eighth of a disc along $\\vec{u}$** as described and motivated just above.\n$$ $$\n - Finally we fit a 1D model with 10 different radiuses along direction $\\vec{u}$:\n\n \\begin{equation} \\tag{5} \\label{lineic_general}\n NO2_{\\vec{u}}(r,t) = f\\big(r, \\vec{v}(t) \\cdot \\vec{u}, T_{\\vec{u}}(r,t), H_{\\vec{u}}(r,t),t\\big)\n \\end{equation}\n$$ $$\n - $\\vec{u}$ [unit vector]: discrete direction that corresponds to the local principal direction of the wind.\n$$ $$\n - $r$ [m]: distance from the cluster barycentre along direction $\\vec{u}$.\n$$ $$\n - $t$ [s]: UTC date of the instantaneous NO2 observation.\n$$ $$\n - $\\vec{v}(t)$ [m.s$^{-1}$]: spatial mean vector, taken close to the barycentre, of the wind at time $t$.\n$$ $$\n - $ NO2_{\\vec{u}}(r,t)$ [mol.m$^{-1}$]: averaged lineic density over directions $\\big[\\vec{u}-\\pi/8,\\vec{u}+\\pi/8\\big]$ of tropospheric NO2 at time $t$ and distance $r$ from the barycentre. \n$$ $$\n - $T_{\\vec{u}}(r,t)$ [°C]: averaged temperature above ground at distance $r$ from the barycentre along $\\vec{u}$ at time $t$.\n$$ $$\n - $H_{\\vec{u}}(r,t)$ [%]: averaged relative humidity above ground at distance $r$ from the barycentre along $\\vec{u}$ at time $t$.\n \n\nIn this way, with further assumptions on the $f$ function, we hope that fitting on our 1D dataset will enable us to distinguish power-plant (and maybe urban activities) contributions to natural (background) contributions expressed by temperature and humidity information (as we mentioned above those two factors are considered as the main natural ones of NO2 emissions)."},{"metadata":{},"cell_type":"markdown","source":"### Let's create the dataframe mentioned above"},{"metadata":{},"cell_type":"markdown","source":"Since wind, temperature and relative humidity information are stored in the same picture and thus have the same coordinates, we can factorize extraction to create the dataframe."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def Create_dataframe(radius_for_surrounding=15, sources_locations=Barycenters):\n '''Whole function to extract data and build the dataframe described above'''\n\n Global_dataframe = pd.DataFrame()\n \n for index, date, hour, wind, temperature, humidity in zip(list(range(len(Time_steps))),\n Time_steps,\n Hours_NO2_instant,\n wind_per_day,\n temperature_per_day,\n humidity_per_day):\n #print('Dealing with index '+str(index)+' over '+str(n_t))\n date_hour = datetime.strptime(date+hour.strftime('/%H/%M/%S'),'%Y/%m/%d/%H/%M/%S')\n \n for location_index in range(len(sources_locations)):\n location = sources_locations[location_index] \n\n \n #let's have local wind direction\n wind_in_surrounding = get_surrounding(location,\n radius_for_surrounding,\n wind, gfs_coords)\n wind_in_surrounding = wind_in_surrounding.reshape(2,-1)\n wind_in_surrounding = wind_in_surrounding[:,np.linalg.norm(wind_in_surrounding,\n axis = 0,\n ord = 2)>pow(10,-2)]\n mean_wind_on_location = np.mean(wind_in_surrounding, axis=1)\n dir_index = vect_to_direction(mean_wind_on_location)\n \n #let's compute lineic functions\n #we extract 8 fields in the same time\n temperature_8_directions = get_8_directions(location,\n radius_for_surrounding,\n temperature,\n gfs_coords)\n humidity_8_directions = get_8_directions(location,\n radius_for_surrounding,\n humidity,\n gfs_coords)\n NO2_8_directions = get_8_directions(location,\n radius_for_surrounding,\n tropospheric_NO2_per_time[index],\n Coords_per_time[index])\n \n #focus on wind direction and its opposite\n for discrete_direction in [dir_index, get_opposite_direction(dir_index)] :\n shapes = []\n \n #since the function get_8_directions computed an average on 2D portions delimited by intervals\n #in the discrete direction (portions of 8th of disk cut in the wind or opposite directions)\n #we will store the middle of each interval in the distance_along_direction array\n distance_along_direction = np.array(temperature_8_directions[discrete_direction][\"upper_bound\"])\n distance_along_direction = distance_along_direction-distance_along_direction[0]/2\n shape_distance_along_direction = distance_along_direction.shape[0]\n shapes.append(shape_distance_along_direction)\n \n #the average temperature over each portion\n density_temperature_1D = np.array(temperature_8_directions[discrete_direction][\"Values\"])\n shape_density_temperature_1D = density_temperature_1D.shape[0]\n shapes.append(shape_density_temperature_1D)\n \n #the average relative humidity over each portion\n density_humidity_1D = np.array(humidity_8_directions[discrete_direction][\"Values\"])\n shape_density_humidity_1D = density_humidity_1D.shape[0]\n shapes.append(shape_density_humidity_1D)\n \n #the average NO2 surfacic density above each portion\n density_NO2_1D = np.array(NO2_8_directions[discrete_direction][\"Values\"])\n shape_density_NO2_1D = density_NO2_1D.shape[0]\n shapes.append(shape_density_NO2_1D)\n \n #since we want equal lengths\n shapes = np.array(shapes)\n shape_0 = np.array([shape_distance_along_direction]*4)\n if(np.not_equal(shapes, shape_0).any()):\n #print(\"Shape issue, continue\")\n continue\n \n #we compute lineic NO2 density from surfacic NO2 density\n #to do so we multiply by the base of the corresponding isosceles triangle\n distance_factor_in_orthogonal_direction = 2*np.mean(distance_along_direction)*np.tan(np.pi/8)\n density_NO2_1D = density_NO2_1D*distance_factor_in_orthogonal_direction\n \n #because dataframe asks for same shape arrays input\n length_arrays = distance_along_direction.shape[0]\n projection_local_wind_along_format = np.array([np.dot(mean_wind_on_location,\n Convention_vectors[discrete_direction])]*length_arrays)\n date_hours_along_format = np.array([date_hour]*length_arrays)\n location_along_format = np.array([location_index]*length_arrays)\n dir_index_along_format = np.array([dir_index]*length_arrays)\n discrete_direction_along_format = np.array([discrete_direction]*length_arrays)\n \n \n #let's create a dataframe\n temp_dict = {\"date\" : date_hours_along_format,\n \"cluster_index\" : location_along_format,\n \"wind_direction_index\" : dir_index_along_format,\n \"discrete_direction\" : discrete_direction_along_format,\n \"distance_along_discrete_direction\" : distance_along_direction,\n \"wind_projection_along_discrete_direction\" : projection_local_wind_along_format,\n \"density_temperature\" : density_temperature_1D,\n \"density_relative_humidity\" : density_humidity_1D,\n \"density_tropospheric_NO2\" : density_NO2_1D}\n\n temp_df = pd.DataFrame(columns=list(temp_dict.keys()), data=temp_dict)\n Global_dataframe = pd.concat([Global_dataframe, temp_df],\n ignore_index=True) \n \n #we convert some data\n Global_dataframe['date'] = pd.to_datetime(Global_dataframe['date'], format='%Y-%m-%d %H:%M:%S')\n Global_dataframe['date'] = np.array(list(map(datetime.timestamp, Global_dataframe['date'])))\n \n dist_in_m = 1e4 #10km=0.086°\n dist_in_lon_lat = 0.0857836897149853\n scale_factor = dist_in_m/dist_in_lon_lat\n \n #from (lon,lat) to meters: distance and NO2 since the latter is now lineic and multiplied by a distance factor\n Global_dataframe['distance_along_discrete_direction'] = scale_factor*Global_dataframe['distance_along_discrete_direction']\n Global_dataframe['density_tropospheric_NO2'] = scale_factor*Global_dataframe['density_tropospheric_NO2'] \n #we export as csv to avoid extracting several times \n print('Saving Dataframe ...')\n Global_dataframe.to_csv(path_or_buf='Global_dataframe.csv',\n index=False)\n return","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Create_dataframe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## V - Fitting: physical analysis to derive an analytical 1D model"},{"metadata":{},"cell_type":"markdown","source":"*Go back to table of contents [⤴](#Content-of-the-rest-of-this-notebook)*"},{"metadata":{},"cell_type":"markdown","source":"Let's extract the dataframe we created:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Global_dataframe = pd.read_csv('Global_dataframe.csv') \nDate_dataframe = pd.read_csv('Date_data.csv')\nTime_steps = Date_dataframe['string_date']\ntimestamp_hours = Date_dataframe['timestamp_date']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### A) Physical discussion for the mathematical form of the surfacic NO2 density"},{"metadata":{},"cell_type":"markdown","source":"#### Physical advection-diffusion paradigm and hypotheses \nWe can assume a certain form for the $f$ function underpinned by physics background, so as to implement and fit additive models:"},{"metadata":{},"cell_type":"markdown","source":"First, for a 2D physical description of total emissions, we write the following formula for a plant isolated (with no big city or airport or any human source of NO2 around):\n\n$$ NO2(\\vec{x},t) = a(t).f_\\text{diff, surf}\\big(\\| \\vec{x}\\|,\\vec{x} \\cdot \\vec{v}(t)\\big) + \\int_{\\vec{x_s}} \\ b\\big(H(\\vec{x_s},t),T(\\vec{x_s},t)\\big).f_\\text{diff, surf}\\big(\\|\\vec{x} - \\vec{x_s}\\|,(\\vec{x}-\\vec{x_s}) \\cdot \\vec{v}(t)\\big) \\text{d}^2\\vec{x_s} $$"},{"metadata":{},"cell_type":"markdown","source":"Where :\n\n- $\\vec{x} \\text{ [m]}$ is the 2D vector that represents the position of the point we consider in regards to the barycentre of the specific cluster of power plants.\n$$ $$\n- $\\vec{x_s} \\text{ [m]}$ is a 2D vector that represents the position of any source of NO2 other than the cluster of power plants (natural background or other activity) in regards to the barycentre of the specific cluster of power plants.\n$$ $$\n- $t \\text{ [s]}$ is the date of the observation : we kept one per day knowing NO2 acquisition time and *GFS* one.\n$$ $$\n- $\\vec{v}(t) \\text{ [m.s}^{-1}\\text{]}$ is the value of the wind above the cluster of plants at time $t$. **We assume that the wind is a constant vector over space around the plant, during the typical lifetime $\\tau$ of NO2 particules**.\n$$ $$\n- $NO2(\\vec{x},t) \\text{ [mol.m}^{-2}\\text{]}$ is the average surfacic density of tropospheric NO2 at position $\\vec{x}$ from the barycentre of the cluster at time $t$.\n$$ $$\n- $a(t) \\text{ [mol]}$ is the amount of molecules of NO2 emitted by the cluster of power-plants that are still 'alive' at time $t$. With equation [(2)](#mjx-eqn-factor) and its related paragraph we discussed the assumptions we have taken on this topic, and how the cluster daily emission can be derived from $a(t)$ knowing $\\tau$.\n$$ $$ \n- $f_{\\text{diff, surf}}\\big(\\| \\vec{x}\\|,\\vec{x} \\cdot \\vec{v}(t)\\big) \\text{ [m}^{-2}\\text{]}$ is the surfacic diffusion function that models the way a source at distance $\\| \\vec{x}\\|$ from the point we consider in a wind $\\vec{v}(t)$ sees its emitted particles moved over space. We define it as a normalized function over space so that $a(t)$ accounts exactly for what we want:\n$$ \\forall \\vec{v}(t) \\int_{\\vec{x_s}} f_\\text{diff, surf}\\big(\\|\\vec{x_s}\\|,\\vec{x_s} \\cdot \\vec{v}(t)\\big) \\text{d}^2\\vec{x_s} := 1 $$\n$$ $$\n- $H(\\vec{x_s},t) \\text{ [%]}$ is the average relative humidity at position $\\vec{x_s}$ at time $t$.\n$$ $$\n- $T(\\vec{x_s},t) \\text{ [°C]}$ is the average temperature at position $\\vec{x_s}$ at time $t$.\n$$ $$\n- $b\\big(H(\\vec{x_s},t),T(\\vec{x_s},t)\\big) \\text{ [mol.m}^{-2}\\text{]}$ is the surfacic density of NO2 sources other than the cluster of plants, we **assume** to depend only on temperature and humidity, not explicitely on time."},{"metadata":{},"cell_type":"markdown","source":"This expression can be thought analogously with heat diffusion in a 2D space and with a specific source highlighted. Though it seems to **rigorously model emissions in an additive way**, it can not be used as such in a computer model. We will make several approximations to get a practical formula. "},{"metadata":{},"cell_type":"markdown","source":"##### How temperature and relative humidity evolve with space ?"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"stats.pearsonr(Global_dataframe[\"distance_along_discrete_direction\"], Global_dataframe[\"density_temperature\"])","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"stats.pearsonr(Global_dataframe[\"distance_along_discrete_direction\"], Global_dataframe[\"density_relative_humidity\"])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As the Pearson correlation coefficients between temperature and position, and humidity and position show, there is no linear correlation between those variables. In addition, we see small local variations in the 2D pictures at the scale we consider. **We will thus assume that $H$ and $T$ do not depend on $\\vec{x_s}$**. Therefore, averaging the other sources term over the integral gives:\n\n$$ \\begin{split}\\int_{\\vec{x_s}} b\\big(H(\\vec{x_s},t),T(\\vec{x_s},t)\\big).f_\\text{diff, surf}\\big(\\|\\vec{x} - \\vec{x_s}\\|,(\\vec{x}-\\vec{x_s}) \\cdot \\vec{v}(t)\\big) \\text{d}^2\\vec{x_s} \\simeq \\\\\n b\\big(H(t),T(t)\\big) \\int_{\\vec{x_s}} f_\\text{diff, surf}\\big(\\|\\vec{x} - \\vec{x_s}\\|,(\\vec{x}-\\vec{x_s}) \\cdot \\vec{v}(t)\\big) \\text{d}^2 \\vec{x_s} \\end{split}$$\n\nwhere $H(t)$ and $T(t)$ are the averages over space of their corresponding smooth function. Then using the normalized property of $f_{\\text{diff, surf}}$ raises: \n\n$$ \\int_{\\vec{x_s}} b\\big(H(\\vec{x_s},t),T(\\vec{x_s},t)\\big).f_\\text{diff, surf}\\big(\\|\\vec{x} - \\vec{x_s}\\|,(\\vec{x}-\\vec{x_s}) \\cdot \\vec{v}(t)\\big) \\text{d}^2\\vec{x_s} \\simeq b\\big(T(t),H(t)\\big)$$\n\nthat leads us to the following simplification for the NO2 model:\n\n$$\\boxed{ NO2(\\vec{x},t) = a(t).f_\\text{diff, surf}\\big(\\| \\vec{x}\\|,\\vec{x} \\cdot \\vec{v}(t)\\big) + b\\big(T(t),H(t)\\big)}$$"},{"metadata":{},"cell_type":"markdown","source":"As a first attempt we assume :\n- $b$ is affine in $T$ and $H$.\n- $f_\\text{diff, surf}$ is a gaussian function distorted by wind."},{"metadata":{},"cell_type":"markdown","source":"### B) Let's go into an analytical 1D model"},{"metadata":{},"cell_type":"markdown","source":"As explained previously, we think that we have better chance to perform a good fit reducing the 2D data in a proper 1D direction. Hence, for a given discrete direction $\\vec{u}$ (that we found to be the closest to the local wind direction), we can reduce our surfacic physical model into a lineic one:\n\n$$NO2_{\\text{lin}}\\big(r,t,\\vec{v}(t)\\cdot\\vec{u},T(t),H(t)\\big) = a(t). f_{\\text{diff, lin}}\\big(r,\\vec{v}(t)\\cdot\\vec{u}\\big) + b_{\\text{lin}}\\big(T(t),H(t)\\big)$$\n$$ \\forall \\vec{v}(t) \\int_{\\mathbb{R}} f_\\text{diff, lin}\\big(r,\\vec{v}(t)\\cdot\\vec{u}\\big) \\text{d}r := 1$$\n\nWith $f_\\text{diff, lin}$ now expressed in m$^{-1}$ and $b_{\\text{lin}}$ in mol.m$^{-1}$. While $a(t)$ stays in mol.\n\nQuick reminder of variables and physical quantities defined in equation [(5)](#mjx-eqn-lineic_general):\n- $\\vec{u}$ [unit vector]: discrete direction that corresponds to the local principal direction of the wind.\n- $r$ [m]: distance from the cluster barycentre along direction $\\vec{u}$.\n- $t$ [s]: UTC date of the instantaneous NO2 observation.\n- $\\vec{v}(t)$ [m.s$^{-1}$]: spatial mean vector, taken close to the barycentre, of the wind at time $t$.\n- $ NO2_{\\text{lin}}\\big(r,\\vec{v}(t)\\cdot\\vec{u}\\big)$ [mol.m$^{-1}$]: modeled lineic NO2 density along direction $\\vec{u}$ of tropospheric NO2 at distance $r$ from the barycentre under a stationnary wind projection $\\vec{v}(t)\\cdot\\vec{u}$. \n- $T(t)$ [°C]: averaged temperature above ground along $\\vec{u}$ at time $t$.\n- $H(t)$ [%]: averaged relative humidity above ground along $\\vec{u}$ at time $t$.\n\nLet's dive into analytical modelisation:\n$$f_{\\text{diff, lin}}\\big(r,\\vec{v}(t)\\cdot\\vec{u}\\big) = \\Big[h\\big[\\vec{v}(t)\\cdot\\vec{u}\\big]*g\\Big](r) $$\n$$b_{\\text{lin}}\\big(T(t),H(t)\\big) = \\alpha T(t) + \\beta H(t) + \\gamma$$\n\nWith $*$ representing convolution over $\\mathbb{R}$, $h$ and $g$ are defined as:\n\n$$h\\big[\\vec{v}(t)\\cdot\\vec{u}\\big](r) = \\frac {\\mathbb{1}_{r\\vec{v}(t)\\cdot\\vec{u} > 0}} {\\tau\\vec{v}(t)\\cdot\\vec{u}} exp\\Big(-\\frac {r} {\\tau\\vec{v}(t)\\cdot\\vec{u}}\\Big) $$ \n$$$$\n$$g(r) = \\frac {1} {\\sqrt{2\\pi}\\sigma}exp\\Big(-\\frac{r^2}{2\\sigma^2}\\Big)$$\n\nSince they are both normalized, by applying Fubini's theorem we derive $f_\\text{diff, lin}$ normalization property. After computing analytically the convolution, we get:\n\n$$f_{\\text{diff, lin}}\\big(r,\\vec{v}(t)\\cdot\\vec{u}\\big) = \\frac{1} {2\\tau\\vec{v}(t)\\cdot\\vec{u}} exp\\Bigg(\\frac {\\sigma^2} {2\\big(\\tau \\vec{v}(t)\\cdot\\vec{u}\\big)^2}\\Bigg) exp\\Bigg(-\\frac {r} {\\tau \\vec{v}(t)\\cdot\\vec{u}}\\Bigg) \\Bigg[ 1+ erf\\bigg(\\frac{1}{\\sqrt{2}\\sigma}\\Big(r-\\frac {\\sigma^2} {\\tau\\vec{v}(t)\\cdot\\vec{u}}\\Big)\\bigg)\\Bigg]$$"},{"metadata":{},"cell_type":"markdown","source":"Thanks to equation [(2)](#mjx-eqn-factor), we know how to derive cluster NOx emissions when obtaining the fit 'parameter' $a$, that accounts for the amount of molecules of NO2 emitted by the cluster of power-plants that are still 'alive' at time $t$. Let's see how to perform such a fit!"},{"metadata":{},"cell_type":"markdown","source":"### C) Fitting implementation"},{"metadata":{},"cell_type":"markdown","source":"We will focus on each cluster separately.\nThe way we want to model $NO2_{\\text{lin, observed}}(r,t)$ is quite particular. Indeed it has several parameters:\n- Constant parameters: $\\tau$, $\\sigma$, $\\alpha$, $\\beta$ and $\\gamma$.\n- **A time-dependent parameter $a(t)$** that we particularly want to know! \n\n\nSince $\\tau$, $\\sigma$, $\\alpha$, $\\beta$ and $\\gamma$ have a physical sense we can provide some rough estimates. Then, for a given set of constant variables $(\\tau,\\sigma, \\alpha, \\beta, \\gamma)$, we are going to compute for each moment $t$ the value $a(t)$ as the following optimum:\n\n\\begin{equation} \\tag{6} \\label{leastsquares} \\boxed{ a(t) = \\text{argmin}_a \\sum_r \\Big[NO2_\\text{lin, mod}[a,\\tau,\\sigma, \\alpha, \\beta, \\gamma]\\big(r,\\vec{v}(t)\\cdot\\vec{u},T(t),H(t)\\big)- NO2_{\\text{lin, observed}}(r,t)\\Big]^2 \\, \\,}\n\\end{equation}\n\nwhere model parameters (written between $[\\ldots]$) are distinguished from variables (written between $(\\ldots)$) as follow:\n\n$$ NO2_\\text{lin, mod}[a,\\tau,\\sigma, \\alpha, \\beta, \\gamma]\\big(r,\\vec{v}(t)\\cdot\\vec{u},T(t),H(t)\\big) = a. f_{\\text{diff, lin}}[\\tau,\\sigma]\\big(r,\\vec{v}(t)\\cdot\\vec{u}\\big) + b_{\\text{lin}}[\\alpha, \\beta, \\gamma]\\big(T(t),H(t)\\big) $$\n\n*NB*: In a next version, we will implement a fit for every parameter, including $\\tau$, $\\sigma$, $\\alpha$, $\\beta$ and $\\gamma$. The reason we \"manually\" set those parameters for now, and focus on $a(t)$ variable is that we can give a physical sense to those former. We will use optimization formalism to add constraints. \n\nIn the following model function, $X$ is the vector of variables we consider :\n$$ $$\n$$X = \\Big(r,t,\\vec{v}(t)\\cdot\\vec{u},T(t),H(t)\\Big)$$\n$$$$\nWe have one vector per sample in the *Global_dataframe* dataset."},{"metadata":{},"cell_type":"markdown","source":"We define the form of the function *model* we want to fit and links between parameters. The function *find_a* does the optimisation process to find $a(t)$ discussed above and returns also the fit relative least-squares error over the 1D segment used. We use this error as rough incertitude approximation."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def f_diff(r, v_dot_u, tau, sigma):\n x0 = tau*v_dot_u\n return (1/(2*np.abs(x0))) * np.exp((1/2)*(sigma/x0)**2) * \\\n np.exp(-r/x0)*(1+erf((1/(np.sqrt(2)*sigma))*(np.sign(x0)*r-sigma**2/x0)))\n\ndef b(T, H, alpha, beta, gamma):\n return alpha*T + beta*H + gamma\n\ndef model(r, t, v_dot_u, T, H, tau, sigma, alpha, beta, gamma, a=2e-2):\n return a*f_diff(r, v_dot_u, tau, sigma) + b(T, H, alpha, beta, gamma)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"def find_a(clust_ind, t, tau, sigma, alpha, beta, gamma, \n data=Global_dataframe, error_threshold=0.2, plot=False):\n '''For given set of variables and a given time t, finds the corresponding\n activity a(t) of the cluster of power-plants, computed as the optimal number that fits :\n N02(x,t) = model(a(t),x,wind,temp,humdity,tau,sigma,alpha,beta,gamma)'''\n \n #we select in our dataframe the cluster, day and we focus on wind direction only\n index = np.logical_and(data[\"cluster_index\"] == clust_ind,\n np.logical_and(data['date'] == t, \n data['wind_projection_along_discrete_direction']>0))#means in wind direction\n \n if(not np.any(index)):#in case no data is consistent for the selection above\n raise ValueError('No data for this specific day, cluster and wind direction')\n return \n \n #values of observed NO2\n expected_values = np.array(data['density_tropospheric_NO2'][index])\n n_dist = len(expected_values)\n \n rs = data['distance_along_discrete_direction'][index]\n v_dot_us = data['wind_projection_along_discrete_direction'][index]\n Ts = data['density_temperature'][index]\n Hs = data['density_relative_humidity'][index]\n \n #diffusion function defined above\n f_diffs = np.array(list(map(lambda rv : f_diff(rv[0], rv[1], tau, sigma),\n zip(rs, v_dot_us))))\n \n #we know that a(t)*f_diff(X) can not be greater than NO2(X,t)\n upper_bound = np.min(expected_values/f_diffs)\n if upper_bound==0: # dealing with bad NaN extrapolation \n raise ValueError('Bad NaN extrapolation')\n return \n \n n_a = 50\n activity_range = np.linspace(0,upper_bound,n_a)#where a can be found\n \n #compute plant contrib and background contrib for each potential a, deduce model value\n models_value = list(map(lambda a : list(map(lambda rvTH : model(rvTH[0], t, rvTH[1],\n rvTH[2], rvTH[3], \n tau, sigma, alpha, \n beta, gamma, a=a),\n zip(rs, v_dot_us, Ts, Hs))),\n activity_range))\n \n errors = list(map(lambda model_values : np.linalg.norm((expected_values-np.array(model_values)), ord = 2),\n models_value))\n \n #select best discrete a to fit expected values with the selected parameters (norm 2 argmin)\n best_index = np.argmin(errors)\n relative_min_error = errors[best_index]/np.linalg.norm(expected_values,ord = 2)\n \n #store the argmin and the corresponding model\n a_opt = np.float(activity_range[best_index])\n model_opt_value = np.array(models_value[best_index])\n \n if(relative_min_error>error_threshold):\n raise ValueError('Relative error too high')\n return \n \n if plot:\n fig = plt.figure()\n plt.xlabel('Distance from the cluster along wind direction (m)')\n plt.ylabel('Lineic density of NO2 (kg/m)')\n plt.title('Daily fitting results for cluster {}'.format(clust_ind))\n plt.scatter(rs*np.sign(v_dot_us), model_opt_value, label='model', marker='x')\n plt.scatter(rs*np.sign(v_dot_us), expected_values, label='observation', marker='+')\n plt.scatter(rs*np.sign(v_dot_us), model_opt_value-gamma, label='cluster-related', marker='+')\n plt.ylim([0, 1.2*np.amax(expected_values)])\n plt.legend()\n plt.show()\n \n return a_opt, relative_min_error","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"# parameters for our model: set manually for every cluster and day \n#(since they are inherent in physics thus not particular)\ntau = 7.2e3#2 hours\nsigma = 2e3\nalpha = 3e-4\nbeta = 1e-4\ngamma = 1e-2\n\nprint(\"Less is more... We fit manually and this is what works well for now:\\n\")\nprint(\"tau: {} seconds\".format(tau))\nprint(\"sigma: {} km\".format(sigma/1e3))\nprint(\"alpha: {} mol.m-1.°C-1\".format(alpha))\nprint(\"beta: {} mol.m-1.%-1\".format(beta))\nprint(\"gamma: {:.2e} mol.m-1\".format(gamma))\n\n#to browse the list of clusters\nCluster_indexes = Global_dataframe[\"cluster_index\"].unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here is an example of fit for a particular date and cluster of power-plants:"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"a_ex, error_ex = find_a(3, timestamp_hours[130], tau, sigma, alpha, beta, gamma, plot = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"print('The fitting attempt you can see above gives us:\\na(t) = {:.4f} mol with a relative error of {:.4f}%'.format(a_ex,100*error_ex))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the next cell, we fit our model per accessible day and per cluster of power-plants to fill a dataframe we will call *Plant_emissions_from_physical_fit* and we will use to compute emissions factors in section VI.\n\nTo select information we believe to be relevant, we only store data when the relative error of the fit is smaller than 20%. This threshold has been set so as to discard aberrant outputs while taking advantage of a maximum of daily values below so as to compute an average emissions factor."},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#we store the relevant data\na_list = []\ndate_list = []\nerror_list = []\nclust_list = []\nprint(\"Fitting per day and cluster and storing data when result is accurate...\")\nfor day in range(len(Time_steps)):\n for cluster_ind in Cluster_indexes:\n try:\n a_,r_error_ = find_a(cluster_ind, timestamp_hours[day], tau, sigma, \n alpha, beta, gamma, error_threshold=0.20, plot=False)\n a_list.append(a_)\n error_list.append(r_error_)\n date_list.append(Time_steps[day])\n clust_list.append(cluster_ind)\n except:\n continue","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"dict_results = {'date': date_list, 'cluster_index': clust_list, \n 'cluster_activity': a_list, 'error': error_list }\nFitting_results = pd.DataFrame.from_dict(dict_results)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#store data in .csv\nFitting_results.to_csv(path_or_buf='Plant_emissions_from_physical_fit.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"for clust_ind in np.sort(Fitting_results['cluster_index'].unique()):\n print(\"{} days kept to fit a(t) for cluster {}\".format(np.sum(Fitting_results['cluster_index']==clust_ind), clust_ind))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## VI - Final results and critics"},{"metadata":{},"cell_type":"markdown","source":"*Go back to table of contents [⤴](#Content-of-the-rest-of-this-notebook)*"},{"metadata":{},"cell_type":"markdown","source":"Now that we have estimated plants-related emissions in the wind direction for every day accessible after NO2 preprocessing and per cluster we are going to compute and display average, montly and marginal emission factors."},{"metadata":{},"cell_type":"markdown","source":"### A) Results and answer to the competition"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Plant_emissions_from_physical_fit = pd.read_csv('Plant_emissions_from_physical_fit.csv') ","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Cluster_inds = np.sort(Plant_emissions_from_physical_fit['cluster_index'].unique())","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#evolution over time per cluster\nMean_per_cluster = []\n\n\nfor clust_ind in Cluster_inds:\n relative_indexes = Plant_emissions_from_physical_fit['cluster_index'] == clust_ind\n cluster_activity = np.array(Plant_emissions_from_physical_fit['cluster_activity'][relative_indexes])\n mean_on_kept_period = (24*3600/tau)*Molar_mass_of_NO2*np.mean(cluster_activity)##mean of daily NO2 mass close to the plant\n Mean_per_cluster.append(mean_on_kept_period)\n cluster_error = np.array(Plant_emissions_from_physical_fit['error'][relative_indexes])*mean_on_kept_period\n plt.errorbar(range(len(cluster_activity)),\n (24*3600/tau)*Molar_mass_of_NO2*cluster_activity,\n cluster_error,\n fmt='x',\n label='daily cluster emissions')\n plt.axhline(y=mean_on_kept_period, label='mean over kept days', color='orange')\n plt.xlabel(\"Days kept after NAN values processing\")\n plt.ylabel(\"Mass of tropospheric NO2: kg\")\n plt.legend()\n plt.title(\"Mass of tropospheric vertical column of NO2 \\n linked to daily emissions, cluster : {}\".format(clust_ind))\n plt.show()\n \nMean_per_cluster = np.array(Mean_per_cluster)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"#distribution per cluster\nfor clust_ind in Plant_emissions_from_physical_fit[\"cluster_index\"].unique():\n relative_indexes = Plant_emissions_from_physical_fit['cluster_index'] == clust_ind\n cluster_activity = np.array(Plant_emissions_from_physical_fit['cluster_activity'][relative_indexes])\n sns.distplot((24*3600/tau)*Molar_mass_of_NO2*cluster_activity)\n plt.xlabel(\"Mass of troposheric NO2 : kg\")\n plt.ylabel(\"Density function\")\n\nplt.title(\"Distribution of mass of tropospheric vertical column of NO2 \\n linked to daily emissions per cluster\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"Clusters = [[16],[2,12,14],[23],[18],[17,3,30],[13],[11],[19],[28,15,10,1]]#one cluster : one list of indexes\n#emissions\nfirst_approx_emission_per_clust = 365*(1+NO2_to_NOx)*Mean_per_cluster\nsns.distplot(first_approx_emission_per_clust/1000)\nplt.title(\"Distribution of NOx emission per cluster\")\nplt.xlabel(\"Mass of NOx : tons\")\nplt.ylabel(\"Density\")\nplt.show()\n\n#electricity generation\ngeneration_per_clust = []\nfor cluster in Clusters:\n generation_inside_clust = np.array(power_plants[\"estimated_generation_gwh\"][cluster])\n generation_per_clust.append(np.sum(generation_inside_clust)) \ngeneration_per_clust = np.array(generation_per_clust)\nprint(\"Annual total production per cluster \\n {} GWh\".format(np.array(generation_per_clust,dtype=np.float32)))\n\nsns.distplot(generation_per_clust)\nplt.title(\"Distribution of electricity generation per cluster\")\nplt.xlabel(\"Electricity generated : GWh\")\nplt.ylabel(\"Density\")\nplt.show()\n\n#emission factor\nfirst_approx_emission_factor_global = first_approx_emission_per_clust/generation_per_clust\nprint(\"Mean per cluster: {:.2f} NOx kg/GWh\".format(np.mean(first_approx_emission_factor_global)))\nprint(\"Documented values are between 100 (300 in the introduction table) \\n and 1000 (4000 in the same table) kg/GWh\")\n#see for other values : https://books.google.fr/books?id=2L-OIrTBiRMC&pg=PA176&lpg=PA176&dq=NO2x+kg/GWh&source=bl&ots=HpOZw0yPip&sig=ACfU3U3i0w4Ya5ex2vfwOwjNMk27JiG7qg&hl=fr&sa=X&ved=2ahUKEwjws4H_vOPnAhVoDWMBHds2BdAQ6AEwAHoECAcQAQ#v=onepage&q=NO2x%20kg%2FGWh&f=false\n\nsns.distplot(first_approx_emission_factor_global)\nplt.title(\"Distribution of emission factor per cluster\")\nplt.xlabel(\"Emission factor : NOx kg/GWh\")\nplt.ylabel(\"Density\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Last answer: average emission factor for the year we consider"},{"metadata":{"trusted":false,"_kg_hide-input":true,"scrolled":true},"cell_type":"code","source":"sum_generation = power_plants[\"estimated_generation_gwh\"].sum()\nprint(\"Mean emission factor for electricity generation in Puerto Rico = {:.0f} NOx kg/GWh\".format(first_approx_emission_per_clust.sum()/sum_generation))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Comments on this new estimation"},{"metadata":{},"cell_type":"markdown","source":"First, this result is about **2.5 times smaller than the one we got in section III**, which was considered as an overestimation for several reasons we discussed above. Here, we believe the refined model has a **deeper physical meaning** and it is bolstered by literature, taking wind, temperature and humidity information into account. Those additional features are some keys to **distinguish power-plant-related emissions from other sources**, taking advantage of the localized aspect of our power-plants sources. \n\n\nNonetheless, we can emphasize some limits:\n- We chose to select emissions downwind, in an eighth of a disc that corresponds to the local and current wind direction. This choice enabled us to fit properly with a 1D model as you can see in section V. Nonetheless, excluding other directions seems questionable when the wind field is not stationary, when its norm is too small or when it tends to be isotropic. Thanks to histograms and density plots we managed to highlight that those particular conditions were not common in our dataset on Puerto Rico.\n- Our fit was bounded by a 15km circle around clusters. This point is not bound to generate downward estimates (since we integrate on all space a function fit on this bounded space), but it represents a limit of information.\n- The constant parameters of our fit were set manually. We believe our rough estimates were reliable, but we will implement an optimization algorithm soon to refine this approach.\n- We still have uncertainties due to the instant to daily emissions computations (see equation [(2)](#mjx-eqn-factor)). Those latter might be much larger than in equation [(3)](#mjx-eqn-errorbar) because of noise amplification with this tiny 1D reduction. However, we think that, deleting the background, we reduced the bias, that lastly entailed overestimation. Hence we might be able to decrease $\\Delta a/a$ through statistical averaging over the number $N_{days, cl}$ of days kept per cluster: \n\n\\begin{equation} \\label{finalerrorbar} \\tag{7}\n \\frac{\\Delta EF_{cl,y}}{EF_{cl,y}} = \\frac{\\Delta e_{cl,y}}{e_{cl,y}} = \\sqrt{\\bigg(\\frac{\\Delta \\tau}{\\tau}\\bigg)^2+ \\bigg(\\frac{\\Delta k_{mass}}{k_{mass}} \\bigg)^2 + \\frac{1}{N_{days, cl}}\\bigg(\\frac{\\Delta a}{a}\\bigg)^2}\n\\end{equation}\n\nwhere $\\big(\\Delta \\tau / \\tau \\big)^2$, $\\big(\\Delta k_{mass} / k_{mass} \\big)^2$ and $\\big(\\Delta a / a \\big)^2$ are averaged over the $Days_{cl}$ sample. We use our previous crude estimations in order to assess our result errobar:\n \n\\begin{array}{|c|c|c|}\n \\hline\n \\Delta \\tau / \\tau \\sim 100\\% & \\Delta k_{mass} / k_{mass} \\sim 20\\% & \\Delta a / a \\sim 150\\% \\\\\n \\hline\n\\end{array}\n$$ $$\n$$ \\boxed{\\frac{\\Delta EF_{cl,y}}{EF_{cl,y}} \\sim \\frac{\\Delta \\tau}{\\tau} \\sim 100 \\% } $$"},{"metadata":{},"cell_type":"markdown","source":"#### Additional answer: monthly average emission factors for clusters having all power plants with month data on emission "},{"metadata":{},"cell_type":"markdown","source":"After browsing power-plants dataset, the only clusters that have all their plants with monthly values for July 2018 to December 2018 are the clusters with plant 16 (cluster 0) alone and plant 23 alone (cluster 2). We have no data for 2019 and for other clusters. We believed extrapolating data from other years or from one plant to another is hazardous. \nAfter going into details for those two particular power plants, we discovered some inconsistent data and had no time to find other data sources. We won't display any result for this section because of this inconsistency."},{"metadata":{},"cell_type":"markdown","source":"### B) Comments and perspectives"},{"metadata":{},"cell_type":"markdown","source":"After *GPPD* correction, NO2 interpolation, data preprocessing and fitting, we managed to get average emissions factors for power-plants in a sub-national region, using remote sensing data and techniques. \n\nOne can notice that the approach is **generic**. Indeed, we did not take into account any external dataset specific to Puerto Rico (except the *GPPD* update for energy generation). Nonetheless, since we added many constraints from physics (which are by essence and assumption universal) in a general mathematical formulation (not inherent in Puerto Rico) we could understandably **distinguish power-plants-related emissions in a way that could be extended to other regions**.\n\nThe crux of the matter was to **consider the amplitude of the distribution of plant emissions a(t) as a function of time and not as a constant** so as to compute a relevant evolution over time that depends on the power plants activity evolution over time and not only on the weather conditions evolution over time. We believe this added difficulties to fit but made more sense to understand NO2 evolution over space and time.\n\nTo apply this framework to another location, one might think about **the way clusters of power-plants are defined**.\nIndeed, this step was crucial and depended a lot on the characteristic distances that can be crossed by NO2 molecules in troposphere before they disapppear.\nBesides, though the specific geography of the island did not have a crucial role in our analysis, a sub-national region with a lot of urban areas and airports might need some model updates. Indeed, the distribution of sources we distinguished from the clusters of power-plants were supposed to be smooth over space. This is one limitation entailed by the assumptions we made on the form of the distribution of the density of NO2. A sum of external located sources could be considered as a new term in the sum of sources, though it would suppose additional complexity and difficulties to train.\n\nOne additional limitation of our submission was the **conbined effect of the space resolution and of the spatial boundaries** on the fit of the model. We could not extract a lot of spatial information for power-plants that were close to the limits of the domain we considered. Thanks to small characteristic distances inherent in NO2 diffusion, we could extract a valuable information anyway. Pushing the boundaries of our dataset could be a way to improve results.\n\nExtra datasets could also be a boon on understanding emission factors. For instance, **airports and cities activities** could be taken into account in a similar way we modeled power-plant activity, although model complexity would increase and physics understanding would be less clear. A compromise between the variety of variables we consider and their degree of freedom in a model has to be taken."},{"metadata":{},"cell_type":"markdown","source":"# References\n

\n[1] US EPA, OAR. 2016. “Basic Information of Air Emissions Factors and Quantification.” Overviews and Factsheets. US EPA. September 2. https://www.epa.gov/air-emissions-factors-and-quantification/basic-information-air-emissions-factors-and-quantification.\n\n[2] Turconi, Roberto, Alessio Boldrin, and Thomas Astrup. 2013. “Life Cycle Assessment (LCA) of Electricity Generation Technologies: Overview, Comparability and Limitations.” Renewable and Sustainable Energy Reviews 28: 555–65. doi:10.1016/j.rser.2013.08.013.\n\n[3] Seinfeld, John H., and Spyros N. Pandis. 2016. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. John Wiley & Sons.\n\n[4] “Mapped: The World’s Coal Power Plants in 2020.” 2020. Carbon Brief. March 26. https://www.carbonbrief.org/mapped-worlds-coal-power-plants.\n\n[5] “What Is the Difference between Electricity Generation Capacity and Electricity Generation? - FAQ - U.S. Energy Information Administration (EIA).” 2020. February 4. https://www.eia.gov/tools/faqs/faq.php?id=101&t=3.\n\n[6] “Coal Plant Capacity Factor United States 2019.” 2020. Statista. Accessed April 5. https://www.statista.com/statistics/744947/capacity-factor-of-coal-power-plants-in-the-us-by-unit-type/.\n\n[7] “Capacity Factor - an Overview.” 2020. ScienceDirect Topics. Accessed April 5. https://www.sciencedirect.com/topics/engineering/capacity-factor.\n \n[8] Julian, Alberto. 2020. “EIA-923 Input, NOx Emissions and EF Reference.” March 20. https://kaggle.com/ajulian/eia-923-input-nox-emissions-and-ef-reference.\n\n[9] Liu, Fei, Steffen Beirle, Qiang Zhang, Steffen Dörner, Kebin He, and Thomas Wagner. 2016. “NOx Lifetimes and Emissions of Cities and Power Plants in Polluted Background Estimated by Satellite Observations.” Atmospheric Chemistry and Physics 16 (8): 5283–98. doi:10.5194/acp-16-5283-2016.\n\n[10] J.H.G.M. van Geffen, H.J. Eskes, K.F. Boersma, J.D. Maasakkers, and J.P. Veefkind. 2019. “TROPOMI ATBD of the Total and Tropospheric NO2 Data Products.” http://www.tropomi.eu/sites/default/files/files/publicS5P-KNMI-L2-0005-RP-ATBD_NO2_data_products-20190206_v140.pdf?fbclid=IwAR2Jg-Eji-AjhBt31vzwBA-96qCD25vUH3Cb64uSfTcKax0KjODAZB_erfg.\n\n[11] Spicer, C. W. 1982. “Nitrogen Oxide Reactions in the Urban Plume of Boston.” Science 215 (4536): 1095–97. doi:10.1126/science.215.4536.1095.\n\n[12] Silvern, Rachel F., Daniel J. Jacob, Loretta J. Mickley, Melissa P. Sulprizio, Katherine R. Travis, Eloise A. Marais, Ronald C. Cohen, et al. 2019. “Using Satellite Observations of Tropospheric NO2 Columns to Infer Long-Term Trends in US NOx Emissions: The Importance of Accounting for the Free Tropospheric NO2 Background.” Atmospheric Chemistry and Physics 19 (13): 8863–78. doi:10.5194/acp-19-8863-2019.\n\n[13] Fischer, E. V., D. J. Jacob, R. M. Yantosca, M. P. Sulprizio, D. B. Millet, J. Mao, F. Paulot, et al. 2014. “Atmospheric Peroxyacetyl Nitrate (PAN): A Global Budget and Source Attribution.” Atmospheric Chemistry and Physics 14 (5): 2679–98. doi:10.5194/acp-14-2679-2014.\n\n[14] Beirle, S., K. F. Boersma, U. Platt, M. G. Lawrence, and T. Wagner. 2011. “Megacity Emissions and Lifetimes of Nitrogen Oxides Probed from Space.” Science 333 (6050): 1737–39. doi:10.1126/science.1207824.\n\n[15] Schindlbacher, Andreas. 2004. “Effects of Soil Moisture and Temperature on NO, NO2 , and N2O Emissions from European Forest Soils.” Journal of Geophysical Research 109 (D17): D17302. doi:10.1029/2004JD004590."},{"metadata":{},"cell_type":"markdown","source":"*Go back to table of contents [⤴](#Content-of-the-rest-of-this-notebook)*"}],"metadata":{"celltoolbar":"Format de la Cellule Texte Brut","cite2c":{"citations":{"7729617/3P8KXVID":{"URL":"https://www.sciencedirect.com/topics/engineering/capacity-factor","accessed":{"day":5,"month":4,"year":2020},"genre":"ScienceDirect Topics","id":"7729617/3P8KXVID","title":"Capacity Factor - an overview","type":"webpage"},"7729617/57US368G":{"author":[{"family":"pablo","given":""}],"id":"7729617/57US368G","title":"putain c'est chian","type":"book"},"7729617/59J7PD5A":{"ISBN":"9781118947401","abstract":"Expanded and updated with new findings and new features New chapter on Global Climate providing a self-contained treatment of climate forcing, feedbacks, and climate sensitivity New chapter on Atmospheric Organic Aerosols and new treatment of the statistical method of Positive Matrix Factorization Updated treatments of physical meteorology, atmospheric nucleation, aerosol-cloud relationships, chemistry of biogenic hydrocarbons Each topic developed from the fundamental science to the point of application to real-world problems New problems at an introductory level to aid in classroom teaching","author":[{"family":"Seinfeld","given":"John H."},{"family":"Pandis","given":"Spyros N."}],"id":"7729617/59J7PD5A","issued":{"day":4,"month":4,"year":2016},"language":"en","note":"Google-Books-ID: n_RmCgAAQBAJ","number-of-pages":"1146","publisher":"John Wiley & Sons","shortTitle":"Atmospheric Chemistry and Physics","title":"Atmospheric Chemistry and Physics: From Air Pollution to Climate Change","title-short":"Atmospheric Chemistry and Physics","type":"book"},"7729617/8XTAVK43":{"URL":"http://www.tropomi.eu/sites/default/files/files/publicS5P-KNMI-L2-0005-RP-ATBD_NO2_data_products-20190206_v140.pdf?fbclid=IwAR2Jg-Eji-AjhBt31vzwBA-96qCD25vUH3Cb64uSfTcKax0KjODAZB_erfg","author":[{"family":"J.H.G.M. van Geffen","given":""},{"family":"H.J. Eskes","given":""},{"family":"K.F. Boersma","given":""},{"family":"J.D. Maasakkers","given":""},{"family":"J.P. Veefkind","given":""}],"id":"7729617/8XTAVK43","issued":{"day":6,"month":2,"year":2019},"title":"TROPOMI ATBD of the total and tropospheric NO2 data products","type":"article"},"7729617/CI6CRDTG":{"URL":"https://www.eia.gov/tools/faqs/faq.php?id=101&t=3","accessed":{"day":5,"month":4,"year":2020},"id":"7729617/CI6CRDTG","issued":{"day":4,"month":2,"year":2020},"title":"What is the difference between electricity generation capacity and electricity generation? - FAQ - U.S. Energy Information Administration (EIA)","type":"webpage"},"7729617/DIEEMF88":{"URL":"https://www.statista.com/statistics/744947/capacity-factor-of-coal-power-plants-in-the-us-by-unit-type/","abstract":"This statistic shows the capacity factor of U.S.","accessed":{"day":5,"month":4,"year":2020},"container-title":"Statista","id":"7729617/DIEEMF88","language":"en","title":"Coal plant capacity factor United States 2019","type":"webpage"},"7729617/EI3Q2XDU":{"DOI":"10.1029/2004JD004590","URL":"http://doi.wiley.com/10.1029/2004JD004590","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Schindlbacher","given":"Andreas"}],"container-title":"Journal of Geophysical Research","container-title-short":"J. Geophys. Res.","id":"7729617/EI3Q2XDU","issue":"D17","issued":{"year":2004},"journalAbbreviation":"J. Geophys. Res.","language":"en","page":"D17302","page-first":"D17302","title":"Effects of soil moisture and temperature on NO, NO2 , and N2O emissions from European forest soils","type":"article-journal","volume":"109"},"7729617/FQ4BNUBG":{"DOI":"10.1126/science.1207824","URL":"https://www.sciencemag.org/lookup/doi/10.1126/science.1207824","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Beirle","given":"S."},{"family":"Boersma","given":"K. F."},{"family":"Platt","given":"U."},{"family":"Lawrence","given":"M. G."},{"family":"Wagner","given":"T."}],"container-title":"Science","container-title-short":"Science","id":"7729617/FQ4BNUBG","issue":"6050","issued":{"day":23,"month":9,"year":2011},"journalAbbreviation":"Science","language":"en","page":"1737-1739","page-first":"1737","title":"Megacity Emissions and Lifetimes of Nitrogen Oxides Probed from Space","type":"article-journal","volume":"333"},"7729617/GK38AZTZ":{"DOI":"10.5194/acp-14-2679-2014","URL":"https://www.atmos-chem-phys.net/14/2679/2014/","abstract":"Abstract. Peroxyacetyl nitrate (PAN) formed in the atmospheric oxidation of non-methane volatile organic compounds (NMVOCs) is the principal tropospheric reservoir for nitrogen oxide radicals (NOx = NO + NO2). PAN enables the transport and release of NOx to the remote troposphere with major implications for the global distributions of ozone and OH, the main tropospheric oxidants. Simulation of PAN is a challenge for global models because of the dependence of PAN on vertical transport as well as complex and uncertain NMVOC sources and chemistry. Here we use an improved representation of NMVOCs in a global 3-D chemical transport model (GEOS-Chem) and show that it can simulate PAN observations from aircraft campaigns worldwide. The immediate carbonyl precursors for PAN formation include acetaldehyde (44% of the global source), methylglyoxal (30%), acetone (7%), and a suite of other isoprene and terpene oxidation products (19%). A diversity of NMVOC emissions is responsible for PAN formation globally including isoprene (37%) and alkanes (14%). Anthropogenic sources are dominant in the extratropical Northern Hemisphere outside the growing season. Open fires appear to play little role except at high northern latitudes in spring, although results are very sensitive to plume chemistry and plume rise. Lightning NOx is the dominant contributor to the observed PAN maximum in the free troposphere over the South Atlantic.","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Fischer","given":"E. V."},{"family":"Jacob","given":"D. J."},{"family":"Yantosca","given":"R. M."},{"family":"Sulprizio","given":"M. P."},{"family":"Millet","given":"D. B."},{"family":"Mao","given":"J."},{"family":"Paulot","given":"F."},{"family":"Singh","given":"H. B."},{"family":"Roiger","given":"A."},{"family":"Ries","given":"L."},{"family":"Talbot","given":"R.W."},{"family":"Dzepina","given":"K."},{"family":"Pandey Deolal","given":"S."}],"container-title":"Atmospheric Chemistry and Physics","container-title-short":"Atmos. Chem. Phys.","id":"7729617/GK38AZTZ","issue":"5","issued":{"day":14,"month":3,"year":2014},"journalAbbreviation":"Atmos. Chem. Phys.","language":"en","page":"2679-2698","page-first":"2679","shortTitle":"Atmospheric peroxyacetyl nitrate (PAN)","title":"Atmospheric peroxyacetyl nitrate (PAN): a global budget and source attribution","title-short":"Atmospheric peroxyacetyl nitrate (PAN)","type":"article-journal","volume":"14"},"7729617/GVBUPCAN":{"DOI":"10.1126/science.215.4536.1095","URL":"https://www.sciencemag.org/lookup/doi/10.1126/science.215.4536.1095","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Spicer","given":"C. W."}],"container-title":"Science","container-title-short":"Science","id":"7729617/GVBUPCAN","issue":"4536","issued":{"day":26,"month":2,"year":1982},"journalAbbreviation":"Science","language":"en","page":"1095-1097","page-first":"1095","title":"Nitrogen Oxide Reactions in the Urban Plume of Boston","type":"article-journal","volume":"215"},"7729617/HDFUDRA3":{"URL":"https://www.epa.gov/air-emissions-factors-and-quantification/basic-information-air-emissions-factors-and-quantification","abstract":"This page includes basic information about air emissions factors and the Compilation of Air Pollutant Emission Factors, AP-42.","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"US EPA","given":"OAR"}],"container-title":"US EPA","genre":"Overviews and Factsheets","id":"7729617/HDFUDRA3","issued":{"day":2,"month":9,"year":2016},"language":"en","title":"Basic Information of Air Emissions Factors and Quantification","type":"webpage"},"7729617/KXG3XUUB":{"URL":"https://kaggle.com/ajulian/eia-923-input-nox-emissions-and-ef-reference","abstract":"Explore and run machine learning code with Kaggle Notebooks | Using data from multiple data sources","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Julian","given":"Alberto"}],"id":"7729617/KXG3XUUB","issued":{"day":20,"month":3,"year":2020},"language":"en","title":"EIA-923 input, NOx emissions and EF reference","type":"webpage"},"7729617/MCJXHUBH":{"DOI":"10.1016/j.rser.2013.08.013","URL":"https://linkinghub.elsevier.com/retrieve/pii/S1364032113005534","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Turconi","given":"Roberto"},{"family":"Boldrin","given":"Alessio"},{"family":"Astrup","given":"Thomas"}],"container-title":"Renewable and Sustainable Energy Reviews","container-title-short":"Renewable and Sustainable Energy Reviews","id":"7729617/MCJXHUBH","issued":{"year":2013},"journalAbbreviation":"Renewable and Sustainable Energy Reviews","language":"en","page":"555-565","page-first":"555","shortTitle":"Life cycle assessment (LCA) of electricity generation technologies","title":"Life cycle assessment (LCA) of electricity generation technologies: Overview, comparability and limitations","title-short":"Life cycle assessment (LCA) of electricity generation technologies","type":"article-journal","volume":"28"},"7729617/NUG9MG77":{"URL":"https://www.carbonbrief.org/mapped-worlds-coal-power-plants","abstract":"Since 2000, the world has doubled its coal-fired power capacity to around 2,045 gigawatts (GW) after explosive growth in China and India. A further 200GW is being built and 300GW is planned.","accessed":{"day":5,"month":4,"year":2020},"container-title":"Carbon Brief","id":"7729617/NUG9MG77","issued":{"day":26,"month":3,"year":2020},"language":"en","shortTitle":"Mapped","title":"Mapped: The world’s coal power plants in 2020","title-short":"Mapped","type":"webpage"},"7729617/VVIITCC9":{"DOI":"10.5194/acp-16-5283-2016","URL":"https://www.atmos-chem-phys.net/16/5283/2016/","abstract":"Abstract. We present a new method to quantify NOx emissions and corresponding atmospheric lifetimes from OMI NO2 observations together with ECMWF wind fields without further model input for sources located in a polluted background. NO2 patterns under calm wind conditions are used as proxy for the spatial patterns of NOx emissions, and the effective atmospheric NOx lifetime is determined from the change of spatial patterns measured at larger wind speeds. Emissions are subsequently derived from the NO2 mass above the background, integrated around the source of interest. Lifetimes and emissions are estimated for 17 power plants and 53 cities located in non-mountainous regions across China and the USA. The derived lifetimes for the ozone season (May–September) are 3.8 ± 1.0 h (mean ± standard deviation) with a range of 1.8 to 7.5 h. The derived NOx emissions show generally good agreement with bottom-up inventories for power plants and cities. Regional inventory shows better agreement with top-down estimates for Chinese cities compared to global inventory, most likely due to different downscaling approaches adopted in the two inventories.","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Liu","given":"Fei"},{"family":"Beirle","given":"Steffen"},{"family":"Zhang","given":"Qiang"},{"family":"Dörner","given":"Steffen"},{"family":"He","given":"Kebin"},{"family":"Wagner","given":"Thomas"}],"container-title":"Atmospheric Chemistry and Physics","container-title-short":"Atmos. Chem. Phys.","id":"7729617/VVIITCC9","issue":"8","issued":{"day":28,"month":4,"year":2016},"journalAbbreviation":"Atmos. Chem. Phys.","language":"en","page":"5283-5298","page-first":"5283","title":"NO x lifetimes and emissions of cities and power plants in polluted background estimated by satellite observations","type":"article-journal","volume":"16"},"7729617/ZF5TZ9BE":{"DOI":"10.5194/acp-19-8863-2019","URL":"https://www.atmos-chem-phys.net/19/8863/2019/","abstract":"Abstract. The National Emission Inventory (NEI) of the US Environmental\nProtection Agency (EPA) reports a steady decrease in US NOx emissions\nover the 2005–2017 period at a rate of 0.1 Tg N a−1 (53 % decrease\nover the period), reflecting sustained efforts to improve air quality.\nTropospheric NO2 columns observed by the satellite-based Ozone\nMonitoring Instrument (OMI) over the US show a steady decrease until 2009\nbut a flattening afterward, which has been attributed to a flattening of\nNOx emissions, contradicting the NEI. We show here that the\nsteady 2005–2017 decrease in NOx emissions reported by the NEI is in\nfact largely consistent with observed network trends of surface NO2 and\nozone concentrations. The OMI NO2 trend is instead similar to that\nobserved for nitrate wet deposition fluxes, which is weaker than that for\nanthropogenic NOx emissions, due to a large and increasing relative\ncontribution of non-anthropogenic background sources of NOx (mainly\nlightning and soils). This is confirmed by contrasting OMI NO2 trends\nin urban winter, where the background is low and OMI NO2 shows a\n2005–2017 decrease consistent with the NEI, and rural summer, where the\nbackground is high and OMI NO2 shows no significant 2005–2017 trend. A\nGEOS-Chem model simulation driven by NEI emission trends for the 2005–2017\nperiod reproduces these different trends, except for the post-2009 flattening\nof OMI NO2, which we attribute to a model underestimate of free\ntropospheric NO2. Better understanding is needed of the factors\ncontrolling free tropospheric NO2 in order to relate satellite\nobservations of tropospheric NO2 columns to the underlying NOx\nemissions and their trends. Focusing on urban winter conditions in the\nsatellite data minimizes the effect of this free tropospheric background.","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Silvern","given":"Rachel F."},{"family":"Jacob","given":"Daniel J."},{"family":"Mickley","given":"Loretta J."},{"family":"Sulprizio","given":"Melissa P."},{"family":"Travis","given":"Katherine R."},{"family":"Marais","given":"Eloise A."},{"family":"Cohen","given":"Ronald C."},{"family":"Laughner","given":"Joshua L."},{"family":"Choi","given":"Sungyeon"},{"family":"Joiner","given":"Joanna"},{"family":"Lamsal","given":"Lok N."}],"container-title":"Atmospheric Chemistry and Physics","container-title-short":"Atmos. Chem. Phys.","id":"7729617/ZF5TZ9BE","issue":"13","issued":{"day":12,"month":7,"year":2019},"journalAbbreviation":"Atmos. Chem. Phys.","language":"en","page":"8863-8878","page-first":"8863","shortTitle":"Using satellite observations of tropospheric NO 2 columns to infer long-term trends in US NO x emissions","title":"Using satellite observations of tropospheric NO 2 columns to infer long-term trends in US NO x emissions: the importance of accounting for the free tropospheric NO 2 background","title-short":"Using satellite observations of tropospheric NO 2 columns to infer long-term trends in US NO x emissions","type":"article-journal","volume":"19"},"undefined":{"DOI":"10.5194/acp-16-5283-2016","URL":"https://www.atmos-chem-phys.net/16/5283/2016/","abstract":"Abstract. We present a new method to quantify NOx emissions and corresponding atmospheric lifetimes from OMI NO2 observations together with ECMWF wind fields without further model input for sources located in a polluted background. NO2 patterns under calm wind conditions are used as proxy for the spatial patterns of NOx emissions, and the effective atmospheric NOx lifetime is determined from the change of spatial patterns measured at larger wind speeds. Emissions are subsequently derived from the NO2 mass above the background, integrated around the source of interest. Lifetimes and emissions are estimated for 17 power plants and 53 cities located in non-mountainous regions across China and the USA. The derived lifetimes for the ozone season (May–September) are 3.8 ± 1.0 h (mean ± standard deviation) with a range of 1.8 to 7.5 h. The derived NOx emissions show generally good agreement with bottom-up inventories for power plants and cities. Regional inventory shows better agreement with top-down estimates for Chinese cities compared to global inventory, most likely due to different downscaling approaches adopted in the two inventories.","accessed":{"day":5,"month":4,"year":2020},"author":[{"family":"Liu","given":"Fei"},{"family":"Beirle","given":"Steffen"},{"family":"Zhang","given":"Qiang"},{"family":"Dörner","given":"Steffen"},{"family":"He","given":"Kebin"},{"family":"Wagner","given":"Thomas"}],"container-title":"Atmospheric Chemistry and Physics","container-title-short":"Atmos. Chem. Phys.","id":"7729617/VVIITCC9","issue":"8","issued":{"day":28,"month":4,"year":2016},"journalAbbreviation":"Atmos. Chem. Phys.","language":"en","page":"5283-5298","page-first":"5283","title":"NO x lifetimes and emissions of cities and power plants in polluted background estimated by satellite observations","type":"article-journal","volume":"16"}}},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.7"}},"nbformat":4,"nbformat_minor":4}