{"cells":[{"metadata":{},"cell_type":"markdown","source":"In this notebook, we deduce that the measured NOx values need to be scaled to account for cloud cover."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport rasterio as rio\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport os\nimport folium\nimport tifffile as tiff\nfrom datetime import datetime, timedelta\nfrom branca.element import Template, MacroElement\nfrom skimage.transform import resize\nimport branca\nfrom tqdm import tqdm\nfrom scipy.stats import pearsonr\nimport geopy.distance","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"s5p_no2_path = '../input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/'\ngfs_path = '../input/ds4g-environmental-insights-explorer/eie_data/gfs/'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"pop_map_file = '/kaggle/input/population/imageToDriveExample.tif'\npop_map = rio.open(pop_map_file).read(1)\np = resize(pop_map, (148,475))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"s5p_bounds = rio.open(s5p_no2_path+os.listdir(s5p_no2_path)[0]).bounds\ngfs_bounds = rio.open(gfs_path+os.listdir(gfs_path)[0]).bounds\ns5p_top, s5p_bot, s5p_lef, s5p_rit = s5p_bounds[1], s5p_bounds[3], s5p_bounds[0], s5p_bounds[2] #will be flipped while reading\ngfs_top, gfs_bot, gfs_lef, gfs_rit = gfs_bounds[3], gfs_bounds[1], gfs_bounds[0], gfs_bounds[2]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gppd = pd.read_csv('../input/ds4g-environmental-insights-explorer/eie_data/gppd/gppd_120_pr.csv')\ngppd = gppd.drop(['system:index', 'capacity_mw', 'country',\n       'country_long', 'estimated_generation_gwh', 'generation_gwh_2013',\n       'generation_gwh_2014', 'generation_gwh_2015', 'generation_gwh_2016',\n       'generation_gwh_2017', 'geolocation_source', 'gppd_idnr',\n       'source', 'url', 'wepp_id', 'year_of_capacity_data'], axis=1)\ngppd['long'] = gppd['.geo'].apply(lambda x: float(x[31:48]))\ngppd['lat'] = gppd['.geo'].apply(lambda x: float(x[50:66]))\ngppd['lat'] = np.where(gppd['lat']<10, gppd['lat']+10, gppd['lat'])\ngppd['Ix']=gppd.index\ngppd['coords'] = list(zip(gppd.lat, gppd.long))\n\nno2_latp = (gppd['lat']-s5p_top)/(s5p_bot-s5p_top)*100\nno2_longp = (gppd['long']-s5p_lef)/(s5p_rit-s5p_lef)*100\nno2_col = np.round(475*no2_longp/100).astype('int')\nno2_row = np.round(148*no2_latp/100).astype('int')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#nox dataframe of dates with corresponding files\nnodf = {}\nfor file in os.listdir(s5p_no2_path):\n    ao = datetime.strptime(file[8:16], '%Y%m%d')\n    nodf[ao] = file","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"nodf = pd.DataFrame.from_dict(nodf, orient='index', columns=['file'])\nnodf['date'] = nodf.index\nnodf = nodf.sort_values(by='date')\nnodf.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#nox, cloud fraction collated by month\nno2 = {}\ncf = {}\nfor mth in range(1,13):\n    no2[mth]=[]\n    cf[mth]=[]\nfor file in nodf.file:\n    ao = datetime.strptime(file[8:16], '%Y%m%d')\n    no2[ao.month].append(np.flipud(rio.open(s5p_no2_path+file).read(1)))\n    cf[ao.month].append(np.flipud(rio.open(s5p_no2_path+file).read(7)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#wind speed, wind angle collated by month\nwind, angle = {}, {}\nfor mth in range(1,13):\n    wind[mth]=[]\n    angle[mth]=[]\n    \nfor dt in tqdm(pd.to_datetime(nodf.date)):\n    file = gfs_path + 'gfs_'+ dt.strftime('%Y%m%d') + '12.tif'\n    gfs_im = rio.open(file)\n    us = gfs_im.read(4)\n    vs = gfs_im.read(5)\n    wind[dt.month].append(np.sqrt(us**2+vs**2))\n    angle[dt.month].append(np.arctan2(vs,us))\n    \nmonth_av_no2 ={}\nmonth_av_cf ={}\nmonth_av_wind = {}\nmonth_av_ang = {}\nfor i in range(1,13):\n    month_av_no2[i] = np.nanmean(np.array(no2[i]), axis=0)\n    month_av_cf[i] = np.nanmean(np.array(cf[i]), axis=0)\n    month_av_wind[i] = np.nanmean(np.array(wind[i]), axis=0)\n    month_av_ang[i] = np.nanmean(np.array(angle[i]), axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"fuels = ['Coal','Gas','Hydro','Oil','Solar','Wind']\ncolors = ['black','gray','green','darkred','orange','blue']\nfudic = {}\nfor k,v in zip(fuels, colors):\n        fudic[k]=v","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let us zoom into the Mayaguez Gas plant, chosen for its remote location next to the sea and relatively low population density of the surroundings which would be expected to have no other sources of NOx. We take a look at the NO2 column and cloud fraction at the plant on various randomly chosen days."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"r = no2_row[16]\nc = no2_col[16]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"all_yr_no2 = np.empty([12,148,475])\nall_yr_cf = np.empty([12,148,475])\nall_yr_ang = np.empty([12,148,475])\nfor i in range(1,13):\n    all_yr_no2[i-1] = month_av_no2[i]\n    all_yr_ang[i-1] = month_av_ang[i]\n    all_yr_cf[i-1] = month_av_cf[i]\nannual_no2_avg = np.nanmean(all_yr_no2, axis=0)\nannual_ang_avg = np.nanmean(all_yr_ang, axis=0)\nannual_cf_avg = np.nanmean(all_yr_cf, axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2, figsize=(20,5))\nnocf = np.multiply(annual_no2_avg, annual_cf_avg)\ndate = '2018-08-08'\nprint (date)\nno2map = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(1)\ncfmap = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(7)\nsns.heatmap(no2map[60:70,180:200] , cmap='Reds', ax=ax[0])\nsns.heatmap(cfmap[60:70,180:200] , cmap='Reds', ax=ax[1])\nax[0].set_title('NO2 Column')\nax[1].set_title('Cloud Fraction')\n#[x.axis('off') for x in ax]\nax[0].grid()\nax[1].grid()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see that the NO2 measured is almost inversely proportional to the cloud fraction when we zoom into a small area of ~5 km x ~10 km."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3, figsize=(20,5))\nnocf = np.multiply(annual_no2_avg, annual_cf_avg)\ndate = '2018-07-01'\nprint (date)\nno2map = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(1)\ncfmap = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(7)\nsns.heatmap(no2map[60:100,:80] , cmap='Reds', ax=ax[0])\nsns.heatmap(cfmap[60:100,:80] , cmap='Reds', ax=ax[1])\nsns.heatmap(p[60:100,:80] , cmap='Blues', vmin=0, vmax=750, alpha=0.2, ax=ax[2])\nplt.scatter(c,r-60, color=fudic[gppd['primary_fuel'].iloc[i]])\na = angle[7][1][r,c]\nfor i in range(3):\n    ax[i].arrow(c, r-60, 5*np.cos(a), 5*np.sin(a), color='black', head_width=2) \nax[0].set_title('NO2 Column')\nax[1].set_title('Cloud Fraction')\nax[2].set_title('Population Density')\n[x.axis('off') for x in ax]\nax[0].grid()\nax[1].grid()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2, figsize=(20,5))\nnocf = np.multiply(annual_no2_avg, annual_cf_avg)\ndate = '2018-10-12'\nprint (date)\nno2map = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(1)\ncfmap = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(7)\nsns.heatmap(no2map[60:100,:80] , cmap='Reds', ax=ax[0])\nsns.heatmap(cfmap[60:100,:80] , cmap='Reds', ax=ax[1])\nplt.scatter(c,r-60, color=fudic[gppd['primary_fuel'].iloc[i]])\na = angle[10][12][r,c]\nfor i in range(2):\n    ax[i].arrow(c, r-60, 5*np.cos(a), 5*np.sin(a), color='black', head_width=2) \nax[0].set_title('NO2 Column')\nax[1].set_title('Cloud Fraction')\n[x.axis('off') for x in ax]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2, figsize=(20,5))\nnocf = np.multiply(annual_no2_avg, annual_cf_avg)\ndate = '2019-02-02'\nprint (date)\nno2map = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(1)\ncfmap = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(7)\nsns.heatmap(no2map[60:100,:80] , cmap='Reds', ax=ax[0])\nsns.heatmap(cfmap[60:100,:80] , cmap='Reds', ax=ax[1])\nplt.scatter(c,r-60, color=fudic[gppd['primary_fuel'].iloc[i]])\na = angle[2][2][r,c]\nfor i in range(2):\n    ax[i].arrow(c, r-60, 5*np.cos(a), 5*np.sin(a), color='black', head_width=2) \nax[0].set_title('NO2 Column')\nax[1].set_title('Cloud Fraction')\n[x.axis('off') for x in ax]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the cloud fraction affects the NO2 measurements inversely. Areas with higher cloud fraction have significantly lower NO2 measurements than adjacent areas with lower cloud fraction, when there are no sources of NO2. Cloud fraction shows a higher than proportional effect at values >0.3."},{"metadata":{},"cell_type":"markdown","source":"Assuming that Measured NOx = Actual NOx*(1- n* CF),\nwe can divide the measured NO2 column by a scaling factor of (1-n* CF) (cloud factor) to get an estimate of the actual NOx."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3, figsize=(20,5))\ndate = '2018-07-02'\nprint (date)\nno2map = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(1)\ncfmap = rio.open(s5p_no2_path+nodf.loc[date]['file']).read(7)\nsns.heatmap(no2map[60:100,:80] , cmap='Reds', ax=ax[0])\nsns.heatmap(cfmap[60:100,:80] , cmap='Reds', ax=ax[1])\nsns.heatmap(np.divide(no2map, 1-0.9*cfmap)[60:100,:80] , cmap='Reds', ax=ax[2])\na = angle[7][2][r,c]\nfor i in range(3):\n    ax[i].scatter(c,r-60, color=fudic[gppd['primary_fuel'].iloc[i]])\n    ax[i].arrow(c, r-60, 5*np.cos(a), 5*np.sin(a), color='black', head_width=2) \nax[0].set_title('NO2 Column')\nax[1].set_title('Cloud Fraction')\nax[2].set_title('Adjusted NO2 Column')\n[x.axis('off') for x in ax]","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}