{row.name}
Fuel: {row.primary_fuel}
Cap: {row.capacity_mw:.0f} MW
Gen: {row.generation_hourly:.0f} MWh/h
\"\n folium.Marker(\n (lat, lon),\n color=colors[row.primary_fuel.split(\",\")[0]],\n icon=icon,\n tooltip=tooltip,\n ).add_to(m)\n return m","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"m = map_plants(plants)\nprint(\"You can mouse over the points to see their name and other details.\")\nm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Exploring Total NO2 Levels\n## Now let's have a look at a year in NO2 levels\nThe NO2 data comes from the [Sentinel-5 Precursor](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S5P_OFFL_L3_NO2) mission, which provides NO2 columns in mol/m2. The readings are provided daily, although there will be some days with too much cloud cover to be useful.\n\nTo make it a bit more understandable, I'm going to convert from mols/m2 to tons. To do this, I'll use the molar mass of NO2 (46.006 g/mol [source](https://pubchem.ncbi.nlm.nih.gov/compound/Nitrogen-dioxide)) and the total surface area of Puerto Rico (9104 km2 [source](https://en.wikipedia.org/wiki/Puerto_Rico)).\n\nAlthough the satellite data includes stratospheric NO2, for this challenge we're mostly interested in *tropospheric* levels of NO2 (the layer closest to the surface of the earth) as these will more directly indicate outputs from power stations and other sources."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def flipud(raster, affine):\n \"\"\"\n Flip raster N-S and fix affine.\n Corrects upside down S5P data.\n Adapted from here: https://github.com/perrygeo/python-rasterstats/issues/98#issuecomment-471136551\n \"\"\"\n \n raster = np.flipud(raster)\n affine = Affine(\n affine.a,\n affine.b,\n affine.c,\n affine.d,\n -1 * affine.e,\n affine.f + (affine.e * (raster.shape[0]-1))\n )\n return raster, affine","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def calc_no2(arr, aff, admin):\n \"\"\"\n Get mean of mol/m2, mutiply by total area of PR.\n Then convert to grams/tons.\n \"\"\"\n \n no2_mol_mass = 46.006 # g/mol source: https://pubchem.ncbi.nlm.nih.gov/compound/Nitrogen-dioxide\n total_area = 9104 * 1e6 # m^2 area of PR source: https://en.wikipedia.org/wiki/Puerto_Rico\n no2_col_mean = zonal_stats(\n admin,\n arr,\n affine=aff,\n nodata=-999,\n )[0][\"mean\"]\n no2_mol = no2_col_mean * total_area\n no2_g = no2_mol * no2_mol_mass\n no2_ton = no2_g / 1e6\n return no2_ton","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def average_s5p_cf(f, band=7):\n \"\"\"\n Get S5P cloud-cover fraction.\n \"\"\"\n \n with rio.open(f) as rd:\n arr = rd.read(band)\n return np.nanmean(arr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def calc_tot_no2(f, admin, band=1):\n \"\"\"\n Calculate total NO2 over an area.\n \"\"\"\n with rio.open(f) as rd:\n arr = rd.read(band)\n aff = rd.transform\n arr, aff = flipud(arr, aff)\n no2_ton = calc_no2(arr, aff, admin)\n return no2_ton","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def s5_to_df(s5p_path, admin, disable_tqdm=False):\n \"\"\"\n Create a DF with time-series S5P data for the main values.\n Total NO2, strat, trop and CF.\n \"\"\"\n \n s5p_data = []\n for f in tqdm(list(s5p_path.iterdir()), disable=disable_tqdm):\n s5p_data.append({\n \"date\": f.stem.split(\"no2_\")[1].split(\"T\")[0],\n \"trop_no2\": calc_tot_no2(f, admin, band=2),\n \"strat_no2\": calc_tot_no2(f, admin, band=3),\n \"cf\": average_s5p_cf(f),\n })\n s5p_df = pd.DataFrame(s5p_data)\n s5p_df[\"date\"] = pd.to_datetime(s5p_df[\"date\"])\n s5p_df = s5p_df.sort_values(by=\"date\")\n s5p_df = s5p_df.set_index(\"date\")\n return s5p_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def plot_total_no2(s5p_df, max_cf=0.3):\n \"\"\"\n Create a plot of total NO2 levels.\n \"\"\"\n \n sns.set_style(\"whitegrid\")\n s5p_df_cloudfree = s5p_df.loc[s5p_df.cf < max_cf]\n no2_strat_ts = s5p_df_cloudfree[[\"strat_no2\"]]\n no2_trop_ts = s5p_df_cloudfree[[\"trop_no2\"]]\n cf_ts = s5p_df[[\"cf\"]]\n\n fig, ax1 = plt.subplots(figsize=(20,10))\n no2_strat_ts.plot(ax=ax1)\n no2_trop_ts.plot(ax=ax1, color=\"green\", linewidth=4)\n ax1.set_xlabel(\"Date\", fontsize=16)\n ax1.set_ylabel(\"Total NO2 (tons)\", color=\"green\", fontsize=16)\n ax1.set_ylim([0, 25])\n ax1.legend([\"Stratospheric\", \"Tropospheric\"], loc=\"upper left\", prop={'size': 16})\n ax1.set_title(\"Annual variation in total NO2 levels\", fontsize=16)\n ax1.xaxis.grid(False)\n ax1.tick_params(axis='both', which='both', labelsize=12)\n\n ax2 = ax1.twinx()\n cf_ts.plot(\n linewidth=0.1,\n color=\"black\",\n ax=ax2,\n legend=True,\n )\n ax2.set_ylabel(\"Cloud cover fraction\", color=\"black\", fontsize=16)\n ax2.set_ylim([0, 1])\n ax2.legend([\"Cloud cover fraction\"], prop={'size': 16})\n\n ax1.xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))\n fig.autofmt_xdate()\n sns.despine(left=True)\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"s5p_df = s5_to_df(s5p_path, admin)\nprint(\"This plot shows troposheric NO2 in green, with days above 30% cloudy excluded.\\nStratospheric NO2 and cloud fraction are included for information.\")\nplot_total_no2(s5p_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Now what?\nThere isn't a huge pattern in tropospheric NO2. So any hopes of correlating with activity levels at a large scale or something else are probably out the window. However, this was an average for the whole of Puerto Rico, so it's probably time to start looking at specific sites."},{"metadata":{},"cell_type":"markdown","source":"# Wind to the Rescue\n\nLet's first have a look at the capital, San Juan, which has two fossil fuel stations and loads of people driving around. Puerto Rico seems super windy, so before going any further with NO2, let's see how that will affect the results.\n\nThe wind data comes from the [Global Forecast System](https://developers.google.com/earth-engine/datasets/catalog/NOAA_GFS0P25), and there are forecasts every six hours for that moment (i.e. the model's forecast of the moment the model runs, which is the most accurate).\n\nThe GFS data provides many values (temperature, humidity, etc) but for now I'm interested in wind: the U component (east-west) and V component (north-south). In order to convert these to magnitude and direction, I'll follow [these](http://colaweb.gmu.edu/dev/clim301/lectures/wind/wind-uv) handy instructions.\n\n*magnitude = sqrt(u2 + v2)* \n*angle = [arctan2](https://en.wikipedia.org/wiki/Atan2)(u, v)*\n\nThen I need to get compass bearing (not mathematical angle) of the direction where the wind *comes from*: \n*bearing = 270 - angle x 180/π* \nand subtract 360 if the result is above 360!"},{"metadata":{"_kg_hide-output":false,"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def rad2com(rad):\n \"\"\"Convert mathematical polar radians to compass direction.\"\"\"\n com = 90 - rad * 180/pi\n com = com if com > 0 else com + 360\n return com","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":false,"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def com2rad(com):\n \"\"\"Convert compass direction to mathematical polar radians.\"\"\"\n rad = (90 - com) * pi / 180\n rad = rad if rad > -pi else rad + 2 * pi\n return rad","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def get_wind_data_point(gfs_path, x=-66.75, y=18, disable_tqdm=False):\n \"\"\"\n Create wind time-series at a point.\n Calcualte magnitude and math + compass direction.\n \"\"\"\n \n gfs_data = []\n for f in tqdm(list(gfs_path.iterdir()), disable=disable_tqdm):\n date = f.stem.split(\"gfs_\")[1]\n with rio.open(f) as rd:\n u, v = next(rd.sample(xy=[(x, y)], indexes=[4, 5])) # 4 and 5 are the u and v bands\n gfs_data.append({\n \"date\": datetime(int(date[:4]), int(date[4:6]), int(date[6:8]), int(date[8:])),\n \"u\": u,\n \"v\": v,\n })\n gfs_df = pd.DataFrame(gfs_data)\n gfs_df = gfs_df.sort_values(by=\"date\")\n gfs_df = gfs_df.set_index(\"date\")\n gfs_df[\"mag\"] = np.sqrt(gfs_df.u ** 2 + gfs_df.v ** 2)\n gfs_df[\"dir\"] = np.arctan2(-gfs_df.v, -gfs_df.u) # negative so we get the FROM direction!\n gfs_df[\"compass\"] = gfs_df[\"dir\"].map(rad2com)\n return gfs_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"x, y = plants.loc[\"San Juan CC\"].xy\ngfs_df = get_wind_data_point(gfs_path, x=x, y=y)\nprint(\"Let's have a quick look at a sample of wind data from this point to make sure our maths made sense!\")\ngfs_df.head(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems to... negative *u* is wind to the west, and negative *v* is wind to the south. So with both negative, that wind is blowing to the south-west, i.e. *from* the north-east. So we expect an angle between 0° (due north) and 90° (due east), and that is just what we see!"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def make_rose(df):\n \"\"\"\n Convert wind data to wind rose data.\n Adapted from https://gist.github.com/phobson/41b41bdd157a2bcf6e14.\n \"\"\"\n \n def speed_labels(bins, units):\n labels = []\n for left, right in zip(bins[:-1], bins[1:]):\n if left == bins[0]:\n labels.append('calm'.format(right))\n elif np.isinf(right):\n labels.append('>{} {}'.format(left, units))\n else:\n labels.append('{} - {} {}'.format(left, right, units))\n return labels\n\n total_count = len(df)\n calm_count = len(df.loc[df.mag == 0.0])\n speed_bins = [-1, 0, 1, 2, 3, 4, 5, np.inf]\n speed_labs = speed_labels(speed_bins, units=\"m/s\")\n dir_bins = np.arange(-7.5, 370, 15)\n dir_labels = (dir_bins[:-1] + dir_bins[1:]) / 2\n \n rose = (\n df.assign(\n speed_bins=lambda df: pd.cut(df[\"mag\"], bins=speed_bins, labels=speed_labs, right=True)\n ).assign(\n dir_bins=lambda df: pd.cut(df['compass'], bins=dir_bins, labels=dir_labels, right=False)\n ).replace({'dir_bins': {360: 0}})\n .groupby(by=['speed_bins', 'dir_bins'])\n .size()\n .unstack(level='speed_bins')\n .fillna(0)\n .assign(calm=lambda df: calm_count / df.shape[0])\n .sort_index(axis=1)\n .applymap(lambda x: x / total_count * 100)\n )\n \n # Fill empty rows in rose DF where no wind was recorded in that direction\n # Otherwise plotting complains\n empty = {k: 0 for k, v in rose.loc[0].to_dict().items()}\n for d in np.arange(0, 360, 15):\n if d not in rose.index:\n s = pd.Series(empty)\n s.name = d\n rose = rose.append(s)\n rose = rose.sort_index()\n \n return rose","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def draw_rose(rose, size=10):\n \"\"\"\n Draw a wind rose.\n Adapted from: https://gist.github.com/phobson/41b41bdd157a2bcf6e14\n \"\"\"\n \n def convert_dir(directions):\n N = directions.shape[0]\n barDir = directions * np.pi/180. - np.pi/N\n barWidth = 2 * np.pi / N\n return barDir, barWidth\n \n palette = sns.color_palette('inferno', n_colors=rose.shape[1])\n bar_dir, bar_width = convert_dir(np.arange(0, 360, 15))\n\n fig, ax = plt.subplots(figsize=(size, size), subplot_kw=dict(polar=True))\n ax.set_theta_direction('clockwise')\n ax.set_theta_zero_location('N')\n\n for n, (c1, c2) in enumerate(zip(rose.columns[:-1], rose.columns[1:])):\n if n == 0:\n # first column only\n ax.bar(bar_dir, rose[c1].values, \n width=bar_width,\n color=palette[0],\n edgecolor='none',\n label=c1,\n linewidth=0)\n\n # all other columns\n ax.bar(bar_dir, rose[c2].values, \n width=bar_width, \n bottom=rose.cumsum(axis=1)[c1].values,\n color=palette[n+1],\n edgecolor='none',\n label=c2,\n linewidth=0)\n\n leg = ax.legend(loc=(0.75, 0.95), ncol=2)\n xtl = ax.set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'])\n ax.yaxis.set_major_formatter(ticker.PercentFormatter(decimals=0))\n ax.tick_params(axis='both', which='both', labelsize=14)\n \n return fig, ax","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To understand the wind data in a more intuitive way, let's create a [wind rose](https://en.wikipedia.org/wiki/Wind_rose). This shows, for each wind direction (remember: where the wind comes *from*!), what percentage of days the wind blew that direction (how far out the petal stretches) and what the wind speed was on those days (the petal's colour)."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"w_mean = gfs_df.mean()[\"mag\"]\nrose = make_rose(gfs_df)\nprint(f\"Mean wind speed is {w_mean:.2f} m/s\")\nfig = draw_rose(rose)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## And what does the wind tell us?\nPuerto Rico is windy! The wind rose shows that the wind in San Juan almost always blows from the east, which confirms our guess that no pollution will be blown into Puerto Rico (assuming the wind is similar all over the island). There wasn't a single day with wind blowing from another direction than easterly, and very few days with less than about 3 m/s of wind (6 knots for the nautically inclined, otherwise 11 km/h).\n\nBased on this, we expect NO2 plumes to be blown out to the west, and in general we don't expect NO2 to hang around for very long before being blown off Puerto Rico altogether.\n\n# Investigating Plumes\n## Let's see what Puerto Rico looks like on more and less windy days\nFinally we're going to look at the raw-ish NO2 satellite data and see if we can intuitively tell what's going on.\n\nThis should tell us something about the plumes of NO2. If we define 'calm' as 3 m/s (not that calm!), we can get a different picture for days above and below that cut-off."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def s5p_arr_timeseries(s5p_path, gfs_df, max_cf=0.3, disable_tqdm=False):\n \"\"\"\n Create a timeseries of S5P arrays showing NO2 column density.\n Output a list of dicts, including wind magnitude so we can filter with that.\n \"\"\"\n \n s5p_arrs = []\n for f in tqdm(list(s5p_path.iterdir()), disable=disable_tqdm):\n date, time = f.stem.split(\"no2_\")[1].split(\"_\")[0].split(\"T\")\n date = datetime(int(date[:4]), int(date[4:6]), int(date[6:8]), int(time[:2]), int(time[2:4]), int(time[4:]))\n mean_wind_mag = gfs_df.loc[gfs_df.index.date == date.date()][\"mag\"].mean()\n with rio.open(f) as rd:\n cf = np.nanmean(rd.read(7)) # cloud\n cf = cf <= max_cf\n arr = rd.read(2) # tropo\n aff = rd.transform\n\n arr = arr * cf # only include relatively cloud free days\n arr, aff = flipud(arr, aff) # fix the upside down data!\n s5p_arrs.append({\"date\": date, \"wind\": mean_wind_mag, \"arr\": arr, \"aff\": aff})\n return s5p_arrs, aff","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":false,"_kg_hide-input":true},"cell_type":"code","source":"print(\"Loading wind data.\")\ns5p_arrs, s5p_aff = s5p_arr_timeseries(s5p_path, gfs_df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"wind_cutoff = 3 # this is pretty arbitrary!\ns5p_wind = [w[\"arr\"] for w in s5p_arrs if w[\"wind\"] >= wind_cutoff]\ns5p_calm = [w[\"arr\"] for w in s5p_arrs if w[\"wind\"] < wind_cutoff]\ns5p_wind_mean = np.nanmean(np.stack(s5p_wind), axis=0)\ns5p_calm_mean = np.nanmean(np.stack(s5p_calm), axis=0)\nprint(f\"Number of windy days: {len(s5p_wind)}\")\nprint(f\"Number of calm days: {len(s5p_calm)}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## And how does it look overlaid on Puerto Rico?\nBefore comparing, let's take a look at the windy mean on a map of Puerto Rico for some context. As expected, big plumes behind San Juan and some of the other power stations, and mostly lower levels elsewhere. Surprisingly little behind the Costa Sur and EcoElectrica stations (south west), which produced the most energy in 2018."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def cm2rgb(hexcm):\n \"\"\"Create an RGB (scale 0-1) version of the provided hex colormap function.\"\"\"\n def rgbcm(x):\n h = hexcm(x).lstrip('#')\n return tuple(int(h[i:i+2], 16)/255 for i in (0, 2, 4))\n return rgbcm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def add_no2_to_map(arr, m):\n \"\"\"\n Add a raster image to a Folium map.\n \"\"\"\n \n arr = arr.copy() * 1e6 # becomes mol/km2 instead of mol/m2\n no2cm = branca.colormap.linear.YlOrRd_09.scale(0, 20)\n no2cm.caption = 'Tropospheric NO2 density in mol/km2'\n no2cm.add_to(m)\n\n folium.raster_layers.ImageOverlay(\n image=arr,\n bounds=[[18.6,-67.3,],[17.9,-65.2]],\n colormap=cm2rgb(no2cm),\n opacity=0.7,\n ).add_to(m)\n return m","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"m = map_plants(plants)\nm = add_no2_to_map(s5p_wind_mean, m)\nm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we know what the different plumes represent, let's compare windy to calm days."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def compare_no2_fields(arr_wind, arr_calm):\n \"\"\"Plot comparison of windy to calm days.\"\"\"\n fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20, 8))\n ax1.imshow(arr_wind, cmap=\"YlOrRd\", vmin=0, vmax=2e-5)\n ax1.set_title(\"Windy\", fontsize=16)\n ax2.imshow(arr_calm, cmap=\"YlOrRd\", vmin=0, vmax=2e-5)\n ax2.set_title(\"Calm\", fontsize=16)\n sns.despine()\n ax1.xaxis.grid(False)\n ax1.yaxis.grid(False)\n ax2.xaxis.grid(False)\n ax2.yaxis.grid(False)\n ax1.axis(\"off\")\n ax2.axis(\"off\")\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"compare_no2_fields(s5p_wind_mean, s5p_calm_mean)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Pollution hangs around when it's less windy!\nUnsurprisingly...The cutoff of 3 m/s is pretty high, so even on the 'calm' days there is quite a plume to the west. Unfortunately there weren't enough real calm days to do otherwise! As will become clear below, this is quite an important part of our analysis. Running this analysis in different areas will require careful consideration of wind distributions. First step: make a wind rose!\n\n## Does this tell us anything about emissions?\n**Yes!**\n\n# Using NO2 Line Density to Determine Emissions\nWorking from the [paper](https://scihub.bban.top/10.1126/science.1207824) *Megacity Emissions and Lifetimes of Nitrogen Oxides Probed from Space* by Beirle et al. (published in Science), I can use the characteristics of the plume on windy days to determine the emissions of NOx in kg/hr! Two useful additional sources: [supporting material](https://science.sciencemag.org/content/sci/suppl/2011/09/21/333.6050.1737.DC1/Beirle.SOM.pdf) and a [follow-up paper](https://www.atmos-chem-phys.net/16/5283/2016/acp-16-5283-2016.pdf).\n\nThey propose to integrate the line-density along a line extending up- and down-wind of the point of interest, taking the average mol/m over a line stretching perpendicular to the wind direction. This can then be fit with the following model function\n\n\n\nwhere\n- *M* is the overall fit\n- *E* is the emissions scaling factor in mol/s - **this is what we want!**\n- *x* is the location of the measurement\n- *X* is the virtual source location (may be up- or down-wind from the actual source)\n- *x0* is the folding distance, which controls how quickly the plume decays\n- *G* is a Gaussian smoothing function and *σ* is it's standard deviation\n\n## So if we can measure this line density and fit the function, we get the emission rate!\nLet's get to it.\n\nI'll go 50 km up- and down-wind, and at each point get the line density for a swath 40 km wide.\nIdeally I should convert to a projected [CRS](https://docs.qgis.org/3.4/en/docs/gentle_gis_introduction/coordinate_reference_systems.html) or use the [Haversine formula](https://en.wikipedia.org/wiki/Haversine_formula) but for now let's just rely on the fact that, near the equator 1 km ~= 0.01 degrees. This will make my life a bit easier, and I'm anyway only sampling at the resolution of the provided raster data, so won't be losing much (or any) accuracy.\n\nFrom inspection of the wind rose from earlier, the median wind direction is directly westerly. Let's quickly check this:","attachments":{"2020-03-20_16-40-35.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAZsAAAB+CAIAAAB9BCsqAAAAA3NCSVQICAjb4U/gAAAgAElEQVR4Xu19PZDqSpam3sbEIqvRxhjQERODvGKtYo0J1DEGmjUGJsZAbaGJNaDXQW2htWC88lBHbARqC71YA9pCbxzUFuqIjUDjDPQ46FnwLPQs9JxBz0K1zt2T4k8S4qcKqm79ZBr3FlL+nPwy8+jkyXNOfvPlyxcCJ4wARgAj8CEQ+E8fohe4ExgBjABGACGAORqeBxgBjMDHQQBztI8zlrgnGAGMAOZoeA5gBDACHwcBzNE+zljinmAEMAKYo+E5gBHACHwcBDBH+zhjiXuCEcAI3ISjue5pIM+9P136q751vybtZ9s+Q92Z118VWNw4RuBFEAhwNNcyumIhnU5nGuNQa64pey+4hjYL8C/b6HYN+yRLc21DUXTr9vRb426DyyCyWF4QREGo8BwDv3nVuUVjrqV31fFNqjpKDkAuC1yBqwii2BAFHv0lcBUFoXUaN8dUz4HqmpqimjftgG1qEs+sIRdFIJirSLp1cviP9h2/wAjsEIC11t2vNVh6coWpqPuJ5RiKrF000cBnwJ9W/VIcWsl1loGnk+Z9DB6napNA7kW/Wm5PQ1VE/px3qrXBIvLVVQ9HtRRB3NV9ZM3bpdroqjrXheedcrX3RIpX036rXi3l87l8sVStt4fz1QlKVvNe+S5xX+37ci0H1buYD+YjuC2H9VJzcqryTbuQsdq6IONyUG8OL6gPqp23sgSRKK9zL/ulBJEoDQLT5USf8SuMwAECq2m7Wuvv1tpiMug18wkilu8FptWiX6/15gelQw+I4O/VsJrPAvO6q/v41GpSLxbz9zCNq0N/9nk7n2+dbWFbYtrM5toX5z5H9+b9op2LAZ8NcLDRcHTZ0jzRxmrSzJWexs+Wo1YpX2oOpst166vFqFfL58udI/xk4fGC4kEjy165OvB1IAK3Ra+U8zPxEx358gV4zn11eI7jAJD54FfsWKWrfjFOxEv9DYlr/hacGMeK4ucYgUMEVsNartQPTdBlJ3/A0dBcLufPfchDHG1UKzY71RQRLw+2Ta9G9WK11y7GYBr7V9qyX0zk2iEhZrWYjiaL9WRfTkdT/2uY+4ndQjjs2HOeLHtA1lZc+LKaTHx8OFgf8PdqexSEbTlqV+v9eUTD83YuUdyu2eD71Xyy66HvzaSVv7/LVYNiKHx7ivd398XOYRtAeZyIRbL41XLDEzfVh3GDr07qvhnq6RLImq+7t1pMJn7ZcDUoJw7yhzt9MUdbDcv+r+dqADJ9PN95ojAbbh7//rQITJv3qerBpiqao32B2ZfKnpaighxt2izWRt5Xd7vtXI1q+epwMSjD+iv6hUBoMhbY7X35gqTFei4ey8H8nvdrpVw23/IvvFE1FQ/Uce0orpcT8PfVcjGf9mvVUzLgvFe6zzV3TG05aubuS9FC7Lx1HzvEbTXpVMvN/mTayQMrDzCpRacEEt0K2NqeqSF2liv34Wk9Vw7vymAYiRCgx8EI4oY6HS8PfdlX82G/Xb4DEXr0BSTFcjGbDUh5y04uPFQHjV3M0bxtPoCzgo/XsFPNpu6LzbMS4EFz+AFGYI2AN7kPhYcjHO3LspcHpcwBA/ShGeBo81YRlgJsK7bbThAIgZ8tv6BpDAKF/0M8LAPrCotoX74gPVy8WO90Imb5op2NoVUXSoNyKnEyFY9JAMMqiAt3+XK5VMznUvGzgsiOqZ1iZ9AHkPxiIZSXw+rd3ZoxATeCzaJf3Ft0iuX1oGyZ2mra8diZ19V5Kx+STYFlHigrl9PhYJdGfiEriNukfheL6ClQFbuvttq9CDHVG6tjIK5H41KO5rHiVL5Wq1WL2VT8DjbV5za04fHGvzECOwTQtA1r59HLYxwNsaKo6b+H9C98xy22YZCFCknO6CRh2Jbr2KLkNlSWmjX0H4m0WEjuM7s2nG+SFBk+rSEZNuOKJtmVqPArgiSpRweVIwLFCl3LPsh7yQNTN34iMg2lK9KQXZfkZPp0MZpXNYLnWIYg6Iam8qhYRLJt+5FMkz4iHU3gVVqeFUizyws6q2q8v3uWSaQFLzuZETVZ5ti0zch6l1sDRmdox4Qu7puzTAv4Qibtr4RKM7Re+K//61/vagNdon2tB3HzqPMTt+5BmmGofzZsRonAgKSoR3sGFPgG0NIk2djh7jrmbEaIormniC40AiPuNWPp2vegZJVkKQO/HI1P/5opOOZYjGh1TRf+FyNwAgGYlo9U0jfZT+T1XlFJmMuIYRybcT6O5hg6UZBhSieBo/1szbSGTjY0liKsrvEDTGOW9rdGUsSjE3FqT9FAnx31Zl36rI3VuS7t3luG7ieLqVQOmehBZUk6QzmGlWTp45kRow6QaUqN7xy6bHUbMsk86EaAE0ETJOnuOuxapuGmmaSlmza3/gZEmIUhjnTAlgA3GKmsKBb8/GzThR1BFJSNwJ2gYdQcxzno8ZEHVIblqH09tjuekRzH7VGhIhCCKWISiUoBsTNIVIFjYt/9UTVsMe1jlkdaxI8xAocIwEJwL560UNz7lkfN/23Ve47mGrrDPqCJCRwtRugPUlobs1A+NI3XJUmapmMGWkABzuCMu2aSIXXDdIVCmPPCeotRmYPVOpZ5yThBI5kWFKlwwIBsXZv5Vhdi3kDaTNNcjtusuDB+7lgqCKZgWMy4wnGSojWYg2q9/tNxYMq7vtnj8Q8xtic3AoKZr/I0Q5mGRWRod6bwwriiahxlyhwvEKoCTG1mWOkC7ScmU2Di3+rmzCEyPgKgHYu4E4NfDigWxC2ZoWO6BcQFkqWqTubeGY9nBHvw9XIB+GSI51A042/InlFdkmXh+3UiuYY2fowXCsw2j2VZj0Q86SGPE0bg6Qgg6Qmm5+UFHds9mMvB0tsNKCjA7rcaN7RXTVU35knrU7liz2dQgMrAiVsCThn3+1f4azFs9yarjfIO9NXDeeA1lAgel67fggLpdNoengYqQ62Abt2vU1xO2sVs2X8e6y8R0p0hnVq2Hjr93ORHXfadengnqruTguWkH+qXBwYcCYzWRwHbRjc6tUmvmGtOArTDj2krG4sB4D4llGeJEjKQ8YqFcEOKtKAebzXttYeLLxMYNjDBWE4GwW6hiiM0bwGSLtGjwaEp2CrulKereQfZceTb09DUCPcV/8YIHEUA8ZqQ3RnKu+jkIqw31nr6qDWyr399MrAYtoqpWDxb602QMnveLtc8frYYdep5ONsiYvflziA0c2Fl3W0tbpGCL5GtrlXhCzgGjaXyrRC7mNTv72rXm4p5/R22SnfI5DeRLVUhwdFA7j4BD+AcOHp5zfvVfDl0sgkGrgfP1siEDolXo2Y2HkvcF0tlZDYbZamALGPj2QObCjj1gKeRFmnIfu0+cZev1ludXq/TLGXvsqWm39p2M0oHuIFtSWrDuL2DmDuA2us1mh1gsRvqJjKzOWtecZajwVQoe1bW8ftSrV6rlotwlF1tDU7aEG96gP/DCBxFABknBcWQ5XQAh+jw8Uzkmz2w7/QXhSl+f9oUM2SPdrTdqBdLMI3bmQksFxszNJQzbFKFHoHtbi3iADSq4jfwbAG2fPUAJwIDEV8PI0gEEbGcy5dbnf5wMhkN+u16KZer9k5LMMiODCTU0XQeNELb1R+FG9j/5rf2v2C44hvyCODn7WKUQVyoA2c5WkSH8SOMwA0QAPOm/IF5U3S9YC2ePWctfg1Hg1bnvdpF1kiLQa3aibAriCb8TTxdTdq1S/yHgsSCWfGg127Wm+3eYGvyekV3juIGFsP1vd/I0RagE9V62BguKvdqOjjtshVVCD/DCNwEAWAj9fM+kot+DXwuo/dgezKu5GhQ0XzY659uZjXp90ZRW7WboPGClcAy7wzmL9jA6arP4AZeVr1oPeC22sWwc2ZoThOA32IEXgmB1bzfCvv0BJpejjrti1Qc30Cxyw8acE6MAEYAI/CWEbhJfLS33EFMG0YAI/CJEMAc7RMNNu4qRuDDI4A52ocfYtxBjMAnQgBztE802LirGIEPjwDmaB9+iHEHMQKfCAHM0T7RYOOuYgQ+PAKYo334IcYdxAh8IgQwR/tEg427ihH48Ahgjka4M1Udn4hm9PbmgK13dYikhhNGACMQRuCzczR0DynXdZ8SRTMM4ev/piCqMMur1uu3jFvECLxxBD43R7MUrqCkFbVCv41hcmab+33TTEXWvaue7bFSYeGO5YKo7TgYyUia6IhsRX9CqLy30UNMBUbghRF4JU/UN9gMXHN1F7tVxLbb9c+LsLmPq/gF7myB0JSHcZgWvWIiHnk93u1owTVhBN4bAp9WRnN0gfu9W1EkJhw8/IU/IeeqpziRTxBjVfM0ZVa38kDKunQYLjvJw/UlpshJ5rka8XuMwCdC4L2x4NvQi27xhYDSBxeaP6l2iNa4veoX7iWOjB3+pPp2mSE6MEFAHO3FoAo3u8+PV4Iyws2i7yvw3PHu4DcYgasRuD4+2tUkfIUKPJZxNIL3JQTB5Qj9VjHhxe9fDJqlfDZXP7iJ9JKKIvN4l0AnUvf56IDe+zIoADq+0jwSQ/zwcyLwGTkaBN4HTdW5u0QumA+jaiKWq7Xb58NvXlBZIAtcHAGcKnRvcVQlcJtJgiB29z1EZcHPMAKfCAH/DcSfZK89Ux7++DNxz3MHt8A9FQC4pI74dux25ejLKm1dko7bjZGMIPHRJDhjuWsn4z/rXc3hjl2ptyaWZPlC4g9/6D5oD3Cl3lM7gPNjBD4aAp+Po40V5XuQagrXMzSCRLcNztAl8ZGJSrMceewlFI7mg3DnZ+WBkMZqN/OrrqLavBCdb9MksDQ2/ofvdFm1udM5I2nEDzECHwyBTySPel2F67FgBOOhm0afg8Jq2mnWSqlYrhO4f+s5Ve3LLPqlXHl9JQrcGHjR3hhdbQj3D+7u0byOAFwaI/CuEfhs1htgFfEjLP8Mm7nyy+SM1XFakCosaRpj1x7DLfJX1ohurxc5OaMonCeVpXn+nvheUcZn6k1mGDjneNzae5zJjV9jBD40ApdzNNc9vWTPvX8TMFqGAQyNoNLp5+uc7C5LUhlhxlQYCvZ8HKlWCpKTzlxn12YbEscUfj8jHMtCngDOTNcRk/xREXhZNU85B9AZGq4Gfhzr13in2rosKRKfAW8F7F713LnqmODzwfJd7HT7XARvUO4iCXM1H3Q6kXeJ74uvpv32uTwXNfaCmZY9MHaAlG3Nr2lltfDf+rtc3nDX+SyyhlU47wRzj+qz7UemzayHCdzvWa52rgLnWT14SiG4s6+Wv7+7u7vP1zrerYnLab9ehCfZUuv0dX++ZuC+6Gq5Wm8269Vq69Aj4ykE7fIuJ70qaApeeP+/GtbuUrn2BzNCRBxmPXhw0R2MbyIei6dy5c5ks7iWw3arPz93WScMxSXWG/NOuXqRLepq0qo2L55Tz5o01xXaLP1YsX8BNNc19Zql5y1YR9cwalDZvS89nGdQmKgOtyjP2/k7uDn+4lGdt3OJrWUyuJOl7uuji8uGR3Y1arV3hAxK8ZfmaOja7+pF936HKX2rv1fTdrW2uU970S/ny+3+cNBv13LwoU4Ud6wHbt2unbI3X3fvLEdbTZpgt37p/cFgi5/dk/DWEFx28t7Kv8W5wFvq2/psILDEn0YecPqvwtGW/Wo2V6p3nn6dO7CkGGJpiA8th7UcsLPLu4yM+GL5/YEO6n62/YQK9k3BYswn4qXB9slrcLTLO/oucoLMmSv1N6LYtF3fMxvPs8c/UEtgd+eMzg/1aI41Hu8tEqyuKFP8WlXtcQOUXNsczzaaHXu2+xPekGyFmTWkc8rsXU2n/3BNVWo0xEqBYTjJ8LRLpipkyG++SXLdmWOpHJ1kGxrEqIDnklDIMA0d1FGZJEXRjKB6sSv8yZpZ3k+Kuk7nFar2q/+kSMrbTP+0VsE9jZ6Z2hBlw3lE/4uiMvbG1TZkUZRkSaxwhcr6mWsZ3UaFzfBdGAMmmWTl2ZmGXJgb5sa2xYFpcqhdojhlrEsFQm8UmALf6I6t8JAda4LmhUL8J1XRbTB2kUjpKeFT3LGq/5RMZ3aq1HQm/fhn9bjp4JqICEwIwEQ1f3JNRRQb6g4P19IbHEPDLCzIO93mwWR2Z7oi8kxG0EyFTydp7v/8X3Ay/uYbGua2C0H7+PQ335Dw1ltptiFkkqxiEd4w8AwHf3oJtS1ISrcrVSA8y2ZRHrR1DMY38Hwmi1pa3NpSIq9m71wMJZJhMzCx94uVKlSSirDt+hHifWx8OWqBTmEwnTSz8WLPY5qwnwG/Qf/XazUfDUBbEEMG7atJp1rMZQMi3Kpfivu2A1d8Jaat/LZmZMgQz3fWkiL6PMfuW1OQHkv7zcISGc/H7qud0Xy5mHRKYKNxV58Emx+WwVcA0sGLK6h8E0X7JY+jgR+Ef6guJy0oWaCgJLuYJGBOAr5e9QkIQ6spcs5Kldqj+bTfap32lFhMBr16Lg62LYsvoBgp5bL50/6n4FfW9rKVm73R4vwe0POWiKXushtjl8t76/mY3ftm9apf9AS+E1Ucw+TLop0D546AjJbI1npT6ADAFduKvpGTeQlAw5wt9ybzSa/ZHi2hlRTyzvMIQa4te8l53irV4PlqMR3UkFXPhnzoCwyOl3/erq3lzMi2TvTtq74aVVPxo0ogWPdxNIH2CWnCUwiIo2m361z0Sqn7GtImIFkvtZbnV71iLErphJx0EqVmu4MmeiiBisPzdgylaSuXOJlgDQWKwFYgkavW6ijVincxED/XXBamCqh94vfZfM3bdGwSUq3s++rF5AlP0u3Cv6sfkBcm9339Hmw42nNZdYCjzdvZwKSZNmHZrT8naOU/Qe+Nvm7xYh0OlZ5wduKp1+Pxu1J7qxQ+NhawGC4y2AuXnyIlXKAbA/jUxfZsKVwAcYtjmERwtF3VqKH1lDwymVHh4PpCs3izmQZxIrszM5y3yrs1hT7MW44GgxNLlXueBnE59/470tZhp97CE7SUU7WQ5LElbFK/zzZD74DnR7GXfV82PgMQtEYwWHXGuGOZFy1R01gk69m2/UimyYM9WpIp0D/JM7rbOHhFkBT5aFuww0gHpMK0aNhi4MnpH5ZpOmnhQa6sZVBJ9mVPi4rQ/VWXbBy3mKAYNk3oFtrz7Ejc25fs9ht+Gi51WvqP//gPwOU09S/0Fj4Jf/mXfxlVOXTyEbYgZyxsokqGn7mmYT5SPo+qNMMmHjVjRqyHgjycDuEqtr/RrsEVTbIrRQIeLgYbKlVRuoZD84pe4ZjTgYUtVTaJBPF9t2s2pEy4rlO/ySSQ44fKdR0XVBHHiTyBSXCWh5oF7YbXzrHJ7M2iAKAZQcjKUteQ2YKljalcEra1liA6qpnho5ZPuiJXutw/pQ31oas0UMypY22dQuTrvbNn9iMVPdAzWbIbWnhoqSQVxV72PfA4mms8NP5EpKum0jApVhmL9IYLAHeCIYnobzJJxxyYBhEJTfjIMhF5TzxCM8EyLYLY7ap9mR3bpVLEHxsNg1M81nuYYJ7EYOoGXu5+OIf5kZHaZU5Lpmn+y7/8S1QNL/7s17/+dT6/Pt0ItbUeiidwm9O0Bj5JXq2X8zFfzRQNk9WOnif7bK411lQFrk5IMhVB0qXt7DtBIkS3q2icahQq6d90FePh2CyIrILOZBKEYe8nges4j7E0Q0fm3j28ApOTkznYKs1XmEYDwCANk1VkusCApi7jzFgxmrxkQTHHjMCL//x3GaOtawJ5auGc7uHXeIu4RcRydPSGQj0oe53aljZvHkYynk0Oj6NZ+vineKErNwoh7pCk6TjMSGgx8P2CmztsmiF0wyQ4JgwDcJsYnT7gQ7YmCl0rnNn3m+KkbiW9e0CnaeJHrWtIzIZlgbZTIzgexDJHe1AZ1WAqmd8IjYopR8ZstGcWyYoB6kjEoJEoEw0JSTMsfYLA7av/7qULMr5elu0Iex28MpEZlol9Z+gmsfGrcB3bSbBc5sn1OuOumWRImCauEJ5Z68pcsys0ujbF8MKD1riAk3nFXFPiZVrWeJp0xIL4azgfkOFjdDl9DFdIdE3QuG8OB2ZjK8Y02IM5u6vxBCZoNZ5aYF4dxyZzVItJXmDFilRxM6JKZ6xKWlaEB1bUovISxMwwKJatdGcFKMY1JKPSPb5wLkfo1XIiz2j4oATbg0gNkg0QpDdj6jjOToJG7CV5yF58FaANqLdh32mpFqO+Z7eIElJHhZUz834bTN2gDDozWE0HQctbpGWNUt1BTMThyTSaBpUt3tEtkcjXe8PJdNRvljY2cctBtbhWKYBZKBwRbE9zkQYCdIzrSuCc4O4wEqKneIEURd9+I/7+/kLqTq9j2falZjbBTiIF4/4AyAtXnioPNlD2S3eb4N+eBj2oqT0O1mLY7k1WG03uaj4czg/zrubT+RNUbF4F807xfnfY/2WFVIAbLd+ufghYV2sOItrb5QClfWIb0Xw5KKXu1nrjQwq3T45jgrBHuK/mcwAfaTR3SHoHWp4H8ZHJ7OnRDmJGIfXj7vQKafDg5MGvr0Zav62melQvbQ9cYAF4S/VIW8f79lXfIL3YjvUgSpajZrHcGozWaQhhCKv1rW0HdO7sweP6ZAAseu5jsVS2WC5Xmz2/RnY1LKd2I4RaT+Xq66mOWJqnvw0cDiA4s4fnAs8CbTXtVXN38VgscZereqdHXxbDZi4B089buAtgW4g/weEbLAuPo2VLZehCrVaues9CyTvkghQ+MQjne2+/t3Z2h4vjgp4sRp06MmWEM7dWf2ehPWqV8sVqvdVu1upr6MEnxMsXz1bbp+3HkLY3ka2uTSZRTMpYKn+5Mf8JktGEyAIJd2XwFYDpsFqM0HkqGtFcbW/U5nEE+GyFdMrBij2fgVqz1ayXy/WT3G9bDkwBwpigV/CFvYsn7ov1wWIxbOXhmwkn7rCGlpNOGc6ziFTRO+I4nMzoCCQLhCLadyKE19hq1KztrLKW/Vpzb/4L7hHNdRsIA5j1zTwy6Ot12uD7sF2Lh22dAPVrv4JFmShvOTbY6Hvj6U/bgxKPUGBB27PdI4TvLWxXy8V8GeBO6yLIiNc7vEdpufBlWS0P8sOJ6VkTuCOUXPs4eNYZWZsnY0A6eawVWfJtP1y015w6LE1/NaqXC5/5RcQ0eXHCwPC1dp2n24uTiBvYILCa1PPbDcEZUJb9UrZ65uB8Hx+NpJJ0iDt6P5OcAmcOiimJoMKCk4Z9HjJ8OmSpD2O+K0O+t5lIpM/44w/EI2iGQqrBt0nwhVTZnm87AepL+sISL5yN8sd+O5gm+8YdQ250jwUtIdmG4lOsPoVkxxwT3Fu5svAphH/GvGRGUgqipGekQrSucAuKDVFN2a5yeIlQALVLIj6C4bJE6qpuVQr0ccghvI6ZluS9JfbxrC/wxpmZYJMOvgymnclEnwZDq2k2E//dDz8jq5QPxdHWxiRwhhc4wHkBmG9cJZUuCCJ7rFLq9AQ/VgyeuzREv3xnWJzozod/RfOSoEHcrIrAHBs1OGbSSLDaOnt89A3IeR8AMNvUTWBSXkpmCplja8FWmF/+9s9IkWYr7Afot9cFp1v4L7/5E7Jct9XoM8WP0lPcD4zAOQQukdHO1fEG3gMXK1xCRpJh74g///DTDAS6E6f1l1R1aR7HMmcEDV5+IDrY5sxNZ85+Zi6tepMPxFP4K8YUIs1YnlgZzo4ReNcIHHqqv+vunCU+wxWQBceaCZzNfW0GsILX1Af+V0xjTICRjQBJMs5aLz2xVWc2Q2F5wcjqmMT+xApxdozA+0Xgg8holw8Asq78/bc/jcEd4SJ72strjshJ0iwodGaKqnZllS10NXGdCSxLG10nk4FrV2ixEYpsElHPyUfgtQTvgVUf22qfLI1fYgQ+FAKfTUYjCFbgwYptphlbvdtLj2eaYajvDZvZX2VnPoBltyhWKmKFkDgJ7Rmfn0zd+AkcmgWefn4duCRG4KMg8Pk4GpERhSzxaGrXROV/0vDT4OkBjhy7MobSJQprATHNMvb5u1FOtWYZ+o9EnBOvFPRONYHfYQTeDwKfkKMRdEUswvW+yjY+3suOlqWqTubegTCam3ZQcISdwRb4tdnm7pj26ZSY6PbRVEXEOrSnY4dLfEQEPiNHIyheatwThnzSc/4mow2hSPVkRRLZ5MwYO46pQ4hYZN5LbUNZoHC6fvntaa26hqL+EC9KD/iU82nA4dwfFoFPydHA1FaEwGt/lm9/8LibKK7GU1Sa0+kKWIlkIK650WAbszRYECK/i110LhSS6VRkrpMTz1Yl1c4+yDw+5DyJE375iRD4IBa2zxgxxxCYwpgbm+GYcs+oK7IIROFx915jKGTK5n4DQ6BFWjcbaShmimmO6FoyE1nFyYeOxqcrVmM8FlFFOGEEMAKAwCeV0aDnFCurD2RXkI75FV47P8BR1ic7gXfj1t2VrXCuMfaOWi1jTEKIv2c0ZeuiaDCKhtnZM8DDRT4uAp9XRvPG1NYqrOg8GBA98DXHGNmjqS6TccdjSthHtruYBHcssZzGarp01A/u4rpwRozAR0Lgk3M0xNR0gdc4XXlHHpGWwgnug4aCoeCEEcAI+BHAHA3PB4wARuDjIPB59WgfZwxxTzACGIEtApij4bmAEcAIfBwEMEf7OGOJe4IRwAhgjobnAEYAI/BxEMAc7eOMJe4JRgAjgDkangNbBGy9UUgnKYpmRf21Qi1h9DECt0UAc7Tb4vl+a7NVuJhbMW1LF4guL2q3DrX7fpHBlL8nBDBHe0+j9YK0uk6aE1i4AoFiRLFAoIsYcMIIvD8EMEd7R2NmG90X2w6S6Ux664LgulSGORrj29Zfjop3NBiY1LeJwCfjaK6lyxWmor7DLRW4oBZEcxtX7QVnk6Pp5MMDe9TDiqJmDZZXrRckAVeNEXguAp+Jo8GlnpoqS3+AkLHvjaNB6CNWsARVPnOh9HOnwb6coz+MOfnUfeQkI2miI7vNqKQAAA/dSURBVLIVfR9o/Pp2cQ0YgZsg8Jk4GlzqyQt85ia4vWolVpfnVVpShRcPhGbpXbsinb2zgBa6ckbnOcV6VRxwYxiBswh8Jo52Foy3mWEm84KRlpTKUc3WE+mGi5DHs418Zc92f8KNoopOVipeRA/bME7GjUvyspQxRU5CN+vhhBF4MwhgjvZmhiKaEEsRGn+mRVmgo98/8alrjXXtocIyguZClDa4Ell88E4bXFMuFOCUM0l6Ka04Z65+pytA0/cPgry9EeaJlODsGIEXQeDL50rLTp6I5XvLd9LrRa8Yh8s424ub0rvo5GKJUrPdmayuqnfRyceIeL5zW+quIgkX/uwIYBntRb4TN6rUlB/++DNcxnmzDeeariRToH8yZjR/ZcjIJC/yiZ//JOGt543GG1dzPQKYo12P4UvV4OhS9wciXuBvHl43maRj6BaqaxPJ8oUE8UP34XWuPr2WXFz+EyCAOdqbHWRHk7WfiBgwjaOmYc+jHUxkbZohTON6rT6wNDZO/KzLKnYxeN5g4FI3RuCzcTT33Zii2ZpqPBJEpnDl1jA8YSxNczmxwWXs8dhyZ7pxFS8igT6CeByr2lXVhInEvzECz0TgM3E0Z6Z3HxSTeDS7iqpv7ReeCdxLF3N0j6GlmAx9q6bGIk3SrELyHPzP8RnzoVAxksfdnS5qN5lh7jBLuwgqnOk1ELjm5hTXJU465Zx7/xr9e69tuDqf/IfvfiZynaVRudmV6Y7tkMntvaG+O5GvQclRC8l/+tNjrNh3NO7G++Nr6MJlPycChzKaPe6KfIHjBVFsiALPwV/wrxzSuYDTdNewT+qWXdtQFN36nLhe2+uZPv4Z6ojT9M3YGVQHdyLveY7vTuRrqKWSHomPY/16tdw1dOCy7xoB8Fbpjh2vC5YOOhGYqRTNVLpbO2/HUGTNOslwNv0PmK8sBrX7RKrUme7tlFajejYWLw/8lkuLfrXcnl5i+DLvVGsDbK50CVTBgWjnvPG5q18E85Prv2WBeeveozXbmt+yWlzXp0FgNW1Xa/01m1gNm9VmfzRfLufAjGKwAiZbHBb9eq13do4Re9iWw9pdLJ5rhdfQalAt+01S5+18/vLJO21mc+2zZHyawbuwo6Na6t1wCbDX9WhNVIcX9g5nwwjsEVgNa7lSf2PzvpoM92bfq14xFiv198LUsl/ON8+Yhe92na7RqPz+B1qQxbQ3P/eJZCWZ2+19HE18cDmeDmRxwT3Q3OxBHfAU9J17pXnefWh89Iioti43Gg2BYzOsoM5ANnYtQy4kv/nmm4w4tp1xg6HSFWUMQT8sQxF5Ns13x6rApikqmeYkYy1t7yF1LWsNIUXdctMZHNdb/YLtQQzV9ZNlhftxqyZwPR8XgZksamlxy2DIDLs/27dth64IvrhWVKGSVIQz8RE2vHLRzsWIWPasuw14EcV8giAqvZgMevVcPJYDb5h5v1bKZfMBOW9UTcWL78bv6Blfz+WgnK+PvC/Jsg9eS6na+sdqWE0RidJgueiVy3tXoXkrC2eY+dZguljOh/VsnDiAZ96GLJBixd51fkrP6M3Ti/RLHkcj7pvzpxfGJT43Ah538IlhOzSWo2a+eLBhXPbyMVhfJzD7izXzd3Rt/EikOdYX3wGcmo2diQNJMyBREGCWOSaSwWgzKEhP2tWUiiZ3icKDaoTECjqTdB+MGcEzgQ+No3JpcRx4FPpRUKzuOzg9s1VJsygKBFEQzRyHJH7U1JnMZAgk3Ja0XwsFmxG66g7ZZDIZIyi+UkB4JiVF0P7b7xTV5oU99u7Onv89HB5uaUSHpzhhBJ6CwMwY2xQKBh9Iti4JDVn//qdHliV0XdyLbRSdTkJcmBnBhLeS2wo2HM2aWY9ELB2wfSLpTMbg07/5I1Xu6TJafrBibdhbkr4Ds009JMNmXIiw2pUOd0kkST06qBwRoJviNZt/SudR3n//93//t3/7t6eWunn+bDb7q1/9alMt2N6TrC5LG4YtK/vmKE6W8unfWLRCH6UhwzLx343RLjOKowUx27bYFbvm0Qpp7kHcBIb89ttvV6vV0ZxPfPE3f/M3f/u3fxtVCAYWbOdcMNfBCSPwFATsmf1IJcPf7WShoRUaXjwYptFQeX3v2AzH9Y82Wi1nOJpXJRiXBasmQeb4mbgTRH7Nz1AeiniMcgikaCDLjnqz7t+NJvsvfvGLv/7rv34KZC+SNx6P7+sFNwTHAmGWOWTm6AtAJlM/yKIiGMfiAVEkGYNx8hPq+2aEPwQoG81wHH20Y9TuvgDir/7qrx4fgdfcJh1X6q05GQQhuk1DuJbPgwBMGfeYaE9mBKHQaEDM6f3X3ptjp76cGxmNLrCJ3307m1lEgd6j6Y7B9S/BIz+XTSJpmo4ZaHcRWIHOuGsmGVKHKIHCgRei4zgxKhMWLAnHkAR5fII2km2o4kbw2baf9tLbGm+A5FHuqlZly7LAb9LMVAogcs0UyRaNLs38XUNQOTCajaDcAYU6zflhR9+WLW+IxIdKs+xFIPzjP/5jRIu3frSlMUJ2v3VTuL4PhkCSTiJVzbFegdhPMQX/XHdsN5ZMRy2kTR0bjgbcQ8qrv5EfdL6LlqKXHLj05+c4W/BzlXSBpbqmBXLCJhP8ZxuqxYgPlKE0dNNNu2OCZX2vYUdLsg97rrgpSKULQphh7etEfyUvWrXBMl/hV0YQc/JvRa5CSAKTdMeK4lS6BSBkJjdMvismqYxSUv9BFLWCujsztuHjQSBMrK4y42U12FVqs69/9DZyW+72yn2DC4krYndsk5mKosLB7bHmd5tNMNk9lgc/xwhEIgDf5tTDDA7JNzscd6apZrLAMWgn6pqKllQUv4yErAAo5qgSDdrYcDTgUBVtTAiVCsMWeDBCSAKf6ioWW1OkgNBFskIFNGAmwaLlOJMyGZmsKJoCR6xJnrXBu4BWVdFPvKmPqQp6H07JjP8cIvz2ut+O/j/Ko/9MO794VjUg2vz93//9xUVpQdUdQZQbnEamC6IsNxjSmakQstoVDdRvOFZJx777DsQ2VZXWrAEGq1IhgAU4JAhxhRArADv8OPGv4DTw9VRTmwuJJcqUuAIvsrZ67JAGaPQ2tjE6uHW+GD+c8RMjALfDZhh97G4izCDrJtiaPWQqApemaK7b3Sm8PJDgZBI4FXsCsAO/TscygWW64IJAp5Fu7DA5hsipBX3NOkHlTya32SI8BV1D4DT+Fe4wCtBpyez//H//+3d/9+WQ+kuegPrpl7/85SU5n5PHBUntn1xlqZ/w14RPxX/95+/BIKI1N0X6Oa1cWcadmRa9vsHT1fjkQ8Y0G0fosBXml7/9M/JvmADZV7aLi386BFyzwcmsevBZjwDC0fiCLujKyRvRTlh2HH8179WawwsiW4NXVbUT9kE4Xu2N3oAdWDbotnWjim9TDbKEJvKdk/B5eSCdNr25iJ7VtFPN3SXi8cRdvtb3+bddVBoyAS2p8vC4XdykDsE3QEQ706VLm8P5Ph8CwE/q550lF/0aOF8en4hr3A491SNY48Ejmpd5QteQcfzx5JqaQTWUymtrw2xNNjnx4IDiOKGv+sZFYckeCbD1A0n4aMtkZq0qsD2rlyuSpVREk2ko3e5DgRr//tcF8cBB4Uzt5y4kJoBGVAXY+mA12hks8etoBGheElxV2XiqR+aBw0eNBJOB/UF+ZDZ4+OE+CNNmrrgz0F8Mm8X7VDwWT2VLrdFJseh1gFjNR4NtGs1PfG+G5QQasisvTRm16oNdr5f9ElzDUgoEHTjX6+WgtvUhPpIV3UUDKf60eo9Uhh9jBK5E4KNxtNWgnK1utkjgep8vtXqD4aBTz6fAyes9ucwv1/vO69y/V4uFn4svwLkKwqiEpsxyPp3OF8BcV4vpZNgfTPdF5oPW2fui1k71eM955ULExW+EwAfjaHDf2i54yKJT94U8mjYh5s1989WVes8eJ+RBe72QFmgemE9iF+UA3oAXbg4UbKlUAlqKJxKJ1H1+r89YjtrtrbJ0MfSFRAjUuRbRbn4B37NhwwU/OQI7641ju9J39dyUFUrQN4o7suC/FS7NZBLEbGe5ujnSDfbOtS2Sqfht6eA9+APYBJz8gj3LDJKdZD2HzBdPFBjlxP70r+BjFWG1/JzWXV3R05K+NYkDcwxGsCuaZbCUO5NYdiyacKq5qdm7kLhhuuLaEofkVIuNbHR9/0qGO26vFlkMP8QIvBACH4ijg54oe9Qtfwm37m5EtNWkmQU/pkQ8FoshyQT+gKO6eOq+GDxwOS3CvDxwKM7ATY47EakQubNYC+jUEhDyYxeME0TYZ20cJ3UI5OaFXcEJI/AWEHivu87VfBrWq8/bufzR2JLzVu5+q19bLeao7KoP8XWRbh7t7w4CsMJt5ol4tu5tu1YQtjLhW/6vNXCgFITjgbsbbJXR8fgmSKhH/KpfjCX8JhmwIX3OztGLXxv372RfCxvcDkYgEoF3ydFWk3Y+QaTKe5EDySC140Zo8065dGDJMqxVvahtq0EpHgr5BnHOSjcRYSIxf8JDz9YrVZs8ochh1sWgWe/Pd89Xy+UKhcPzKxWRiPaMZjbkrcPBHTaMn2AEXh2B59mjvdAO+KJqXbMrj5M8k/jxD5XK7pZIR5PGbLQRmjOWZVfshi1ZZmM34zmTgZsWGfSGJVyjqxGcuI8DB5cofRWPpExDKiV+VGXj2VZpji4KWrKQBDM4SIauSqIyI5Ay0B5vo99ZXVF2eOXQ9/b0gLiGov4QL0oPTJRryemy+C1G4IUQeHUeem2DYJHgbRUHZdDgxPOdOapw2vIZoflbWE3gxNNniLZcbmwT5u1y0xN9QDw5iCF7IxHm2q565WEzHXuu3clyVIfLJ4JpI46upr1yvlhrdzrNaqnaGj1dEQbnyvHD3fpN+owrwQg8F4F3uetcdxb2mbAlgzU1XYLb01ZJ5gcC7pgplZr94Widhn24cgZ02BAEcdEv59aHCCswO4X913I08Lt1AZtLFLe3Ocw7+URwi/tctJ9XDnha/KXUeKsTRr4nqfX25YD9yUz4JUbgtRF4xxwNNGATUP8QsdT9/eENViDcdIqe2b0vQUTz5ReI8x8j4rmtrdq0lQWXx1zVr5RDhwHXijC3HMlFv5RKlfyq/VvW/oy6FiAhvxSTfQY5uAhGYIvAQeyNF9rcvlS1M4VlfmuyPUvjLzcSQ7HFk77QNy5oyfx38+6JfTPXwrumVOBURh1LX19p5Y4lltNYTZeiwva+1EjjejEClyDw3jkaxJtURZWWxK+/0C+B+5o8EFKNf3Bl7atEF/IRbimc4D5ovussrukVLosRuCkC75+j3RQOXBlGACPwrhF4f9Yb7xpuTDxGACPwoghgjvai8OLKMQIYgVdFAHO0V4UbN4YRwAi8KAKYo70ovLhyjABG4FUR+P/w3C4WJIg9XAAAAABJRU5ErkJggg=="}}},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"median_wind_direction = gfs_df.loc[gfs_df.mag > 3].compass.median()\nprint(f\"Median wind direction: {median_wind_direction:.0f}°\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So this is actually slightly south of exactly west, but let's stick to a east-west cross-section for our analysis. I suspect if I took into account the higher impact of windier days (which may have a stronger than linear impact) I'd end up concluding the median wind *effect* is due west. Why? Because the NO2 plume is going straight east! Note: I could also rotate the NO2 fields on each day to line up with the wind direction for that day, but as the wind is already basically directly east-west, this probably isn't necessary. In a location with more varied wind directions, one could either calculate a separate profile for each direction, or do this NO2 field rotation.\n\nI want to quickly codify this in our list of power plants. However, to avoid double work, I'm going to take another step in simplifying our plants list by merging those that are right next to each other, as they'll be indistinguishable in terms of satellite imagery.\n\nThe table now has five entries, and those where multiple plants were close together have both names listed as well as both fuel types. I've also entered the primary wind direction for our first investigation."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def merge_close_plants(plants, dist_deg=0.05):\n \"\"\"\n Merge plants that are close together.\n Sum numerical values, list others.\n \"\"\"\n \n buff = plants.copy()\n buff.geometry = buff.geometry.buffer(dist_deg)\n buff = gpd.GeoDataFrame(geometry=[p for p in buff.unary_union])\n buff[\"centroid\"] = buff.geometry.centroid\n\n plants = gpd.sjoin(plants, buff, op=\"within\")\n\n grouped = plants.groupby(\"index_right\").sum()\n grouped[\"geometry\"] = plants.groupby([\"index_right\"])[\"centroid\"].first()\n grouped[\"primary_fuel\"] = plants.groupby([\"index_right\"])[\"primary_fuel\"].apply(lambda x: \", \".join(x))\n grouped[\"name\"] = plants.reset_index().groupby([\"index_right\"])[\"index\"].apply(lambda x: \", \".join(x))\n grouped = grouped.set_index(\"name\")\n grouped[\"capacity_factor\"] = grouped[\"generation_hourly\"] / grouped[\"capacity_mw\"]\n \n grouped = gpd.GeoDataFrame(grouped)\n grouped[\"xy\"] = grouped.geometry.apply(lambda geom: [float(g[0]) for g in geom.coords.xy])\n \n return grouped","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plants = merge_close_plants(plants)\nplants = pretty_cols(plants)\nplants[\"wind_dir\"] = \"\"\nplants.at[\"San Juan CC, Palo Seco\", \"wind_dir\"] = 90\nplants[[\"Primary Fuel\", \"Capacity (MW)\", \"Hourly Generation (MWh)\", \"Capacity Factor\", \"wind_dir\"]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we're finally ready to start our in-depth investigation!"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"def coords2pixel(aff, x, y):\n \"\"\"\n Convert x, y in degrees to pixel coordinates in the array.\n Returned in row, col format as expected by numpy.\n Useful reference: https://www.perrygeo.com/python-affine-transforms.html\n \"\"\"\n \n col, row = ~aff * (x, y)\n col, row = int(col), int(row)\n return row, col","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"def get_line_density(arr, aff, x, y, length=100, width=20, angle=0):\n \"\"\"\n Calculate line density plot for given image and x, y coordinates.\n Length and width are in km, angle is mathematical angle.\n (i.e. N=pi/2, E=0, S=-pi/2, W=pi=-pi)\n \"\"\"\n \n arr_lines = np.zeros_like(arr)\n row1, col1 = coords2pixel(aff, x, y)\n samples = []\n up = int(length*1e3 / 2)\n for x_move in range(-up, up + 1, 100):\n line = []\n side = int(width*1e3 / 2)\n for y_move in range(-side, side + 1, 100):\n x_new = x + (x_move * cos(angle) - y_move * sin(angle)) * 1e-5\n y_new = y + (y_move * cos(angle) + x_move * sin(angle)) * 1e-5\n row, col = coords2pixel(aff, x_new, y_new)\n if row < 0 or col < 0:\n continue\n try:\n line.append(arr[row, col])\n arr_lines[row, col] = 1\n arr_lines[row, col] = 1\n except IndexError:\n pass\n if len(line) > 0:\n line_mean = np.mean(line) * width * 1e3 # to go from mol/m2 to mol/m\n samples.append({\"x\": x_move, \"val\": line_mean})\n line_df = pd.DataFrame(samples)\n line_df = line_df.set_index(\"x\")\n return line_df, arr_lines","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":false,"_kg_hide-input":true},"cell_type":"code","source":"x, y = plants.at[\"San Juan CC, Palo Seco\", \"xy\"]\nline_df, sample_loc = get_line_density(s5p_wind_mean, s5p_aff, x, y, length=160, width=40, angle=0)\nxx = np.array(line_df.index.to_list())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The plot below shows the extents of sampling, up- and down-wind and to either side."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def draw_sample_area(arr_no2, arr_sample):\n \"\"\"Plot NO2 data with sample area overlaid.\"\"\"\n fig, ax = plt.subplots(figsize=(12, 6))\n ax.imshow(arr_no2, cmap=\"YlOrRd\", vmin=0, vmax=2e-5)\n ax.imshow(arr_sample, cmap=ListedColormap([(0, 0, 0, 0), (0, 0, 0, .2)]))\n ax.grid(False)\n\n ax.set_title(\"NO2 with sample area\", fontsize=16)\n sns.despine()\n ax.xaxis.grid(False)\n ax.yaxis.grid(False)\n ax.axis(\"off\")\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"draw_sample_area(s5p_wind_mean, sample_loc)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## And this is what the profile looks like!"},{"metadata":{"_kg_hide-output":true,"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_density(measured, fit_xx=None, fit_yy=None, size=12):\n \"\"\"Plot NO2 density versus profile position.\"\"\"\n \n sns.set_style(\"whitegrid\")\n fig, ax = plt.subplots(figsize=(size, size/2))\n x = measured.index\n spl = UnivariateSpline(x, measured.val)\n spl.set_smoothing_factor(0.005)\n ax.plot(x, spl(x), label=\"Measured\", linewidth=4)\n ax.set_ylabel(\"NO2 line density (mol/m)\", fontsize=16)\n ax.set_xlabel(\"x (km)\", fontsize=16)\n ax.set_title(\"Line density profile (positive x is upwind)\", fontsize=16)\n ax.set_xlim([-80e3 ,80e3])\n ticks = ticker.FuncFormatter(lambda x, pos: '{0:g}'.format(x*1e-3))\n ax.xaxis.set_major_formatter(ticks)\n ax.set_ylim([0, 0.5])\n ax.arrow(20e3, 0.1, -40e3, 0, head_width=0.03, head_length=5e3, fc='k', ec='k')\n ax.text(-18e3, 0.12, \"wind direction\", fontsize=20)\n \n if fit_xx is not None:\n ax.plot(fit_xx, fit_yy, c=\"red\", linewidth=4)\n ax.legend([\"Measured\", \"Fit\"])\n \n ax.xaxis.grid(False)\n ax.tick_params(axis='both', which='both', labelsize=12)\n sns.despine(left=True)\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot_density(line_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## So let's set up the model function\nThese are the functions *e*, *G* and *M* explained before."},{"metadata":{"trusted":true},"cell_type":"code","source":"e = np.vectorize(lambda x, x0, X: np.exp((x+X)/x0) if x <= -X else 0)\nG = lambda x, std: np.exp(-(x**2)/(2*(std**2)))/(np.sqrt(2*np.pi)*std)\nM = lambda x, E, B, X, x0, std: E * np.convolve(e(x, x0, X), G(x, std), mode=\"same\") + B","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Now we use a curve fit to find the parameters\nThis will give the values of *E*, *B*, *X*, *x0*, and *σ* that make the model best fit the data."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"def fit_model(M, xx, yy):\n p, cov = curve_fit(\n f=M,\n xdata=xx,\n ydata=yy,\n p0=(30, 0.25, 1e3, 5e4, 1e4),\n bounds=(\n (0, 0, -np.inf, -np.inf, -np.inf),\n (np.inf, np.inf, np.inf, np.inf, np.inf)\n ),\n method=\"trf\",\n )\n return p","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-output":false,"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"E, B, X, x0, std = fit_model(M, xx, line_df.val.to_list())","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\n f\"And these are our values!\\nE:\\t{E:.2f} mol/s\"\n f\"\\nB:\\t{B:.3f} mol/m\"\n f\"\\nX:\\t{X:.0f} m\"\n f\"\\nx0:\\t{x0:.0f} m\"\n f\"\\nstd:\\t{std:.0f}\"\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's plot this model fit against the measured data and see how it looks."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"yy = M(x=xx, E=E, B=B, X=X, x0=x0, std=std)\nplot_density(line_df, xx, yy)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## It looks pretty good!\nNow let's convert this *E* value to something more intelligible than mol/s and check some other things.\n\nFirstly, the NO2 decay time, *τ*, which is by dividing the *x0* distance by the mean wind speed."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"sec_per_hour = 3600\nno2_mol_mass = 46e-3 # kg/mol\n\ntau = (x0 / w_mean) / sec_per_hour # to convert to hours\nprint(f\"NO2 decay time: {tau:.2f} hours\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"which agrees pretty well with the literature range of 3 - 8 hours.\n\nThen we convert *E* in mol/s to output in kg/hr"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"Ekghr = E * sec_per_hour * no2_mol_mass\nprint(f\"NO2 output is {Ekghr:.0f} kg/hr\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The total hourly electricity generation of the San Juan area is 416 MWh/h. We can use this with the emissions above to calculate the emissions factor in kg/MWh!"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"emission_fac = Ekghr / plants.loc[\"San Juan CC, Palo Seco\", \"generation_hourly\"]\nprint(f\"Naive emissions factor: {emission_fac:.2f} kg/MWh\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Is that it!??\nThis is quite above the expected range of 2 - 3 kg/MWh, but I've forgotten one thing...\n\n## Cars!!\nIf I can find a way to account for cars (and other sources of NO2), I should be able to calculate the real value! According to [this random source](http://www.apis.ac.uk/overview/pollutants/overview_nox.htm) from the UK, road traffic can account for more than 30% of emissions.\n\nBut how do we calculate which portion comes from cars and which from other sources?\n\n# Modelling All the Emissions Sites\n## The plan\nI'll select a bunch of sites and then perform the curve-fitting procedure on each one. Except this time, I'll add an extra factor to account for other sources of emission. I'm also going to include a site with *no* power generation, so help solve for the unknown. I'll take total [night time lights](https://www.nasa.gov/topics/earth/earthday/gall_earth_night.html) emissions in a radius around each point and take that as an indicator of non-power-plant NO2 emissions. Then I can disaggregate our *E* factor into two terms:\n\n*E = E_elec + f_other x NTL*\n\nwhere\n- *E* is the total emissions\n- *E_elec* is the emissions just from the power plants (what we're after!)\n- *f_other* is some multiplier of emissions per light-output that we need to discover\n- *NTL* is the known total light emissions for the area"},{"metadata":{},"cell_type":"markdown","source":"## Let's choose our points and see what happens!\nI selected the following points. These are the main fossil-fuel power plants, plus a large city near San Juan called Caguas. For this city, we know that *E_elec* will be exactly zero.\nUsing these locations, along with some knowledge of the fuels and generation activity of each plant, we have enough data points to solve the different unknown *E_elec* values along with the constant *f_other* factor that accounts for non-power NO2 emissions.\n\nI'm also going to exclude the Cambalache and Mayagüez plants: they are both very small producers (less than 1.5% of our selected sites) and the first is directly under San Juan's plume and won't be detectable, and the second is hovering on the edge of our satellite image, so we won't see its full plume. (It would be trivial to get more data to deal with this, but we're trying to keep things as simple as possible for now!)\n\nSo let's see the table with the latest chops and changes."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plants = plants.loc[plants.generation_hourly > 100]\nplants.loc[\"Caguas\"] = {c: 0 for c in plants.columns}\nplants.at[\"Caguas\", \"xy\"] = [-66.04, 18.24]\nplants.at[\"Caguas\", \"geometry\"] = Point([-66.04, 18.24])\nplants.at[\"Caguas\", \"Primary Fuel\"] = plants.at[\"Caguas\", \"primary_fuel\"] = \"none\"\nplants.at[\"Caguas\", \"wind_dir\"] = \"\"\nplants[[\"Primary Fuel\", \"Capacity (MW)\", \"Hourly Generation (MWh)\", \"Capacity Factor\", \"wind_dir\"]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(f\"Still our included plants account for {(100 * plants.capacity_mw.sum() / gppd.capacity_mw.sum()):.0f}% of the total system capacity.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Great! So let's stick a 5 km buffer around each of our chosen sites and see what it looks like (in red below). Note also the new green marker at the town of Caguas."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plants_buff = plants.copy()\nplants_buff.geometry = plants_buff.geometry.buffer(0.05)\n\nm = map_plants(plants)\nfolium.features.GeoJson(\n plants_buff.to_json(),\n style_function=lambda x: {\n \"fillOpacity\": 0.5,\n \"weight\": 0,\n \"fillColor\": \"red\",\n },\n).add_to(m)\nm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Now we get the total night lights in each circle\nFirst our night time lights. Seems to make sense - hot spots around San Juan and other towns."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"with rio.open(viirs_path) as rd:\n viirs = rd.read(1)\n viirs[np.isnan(viirs)] = 0\n viirs_aff = rd.transform\nfig, ax = plt.subplots(figsize=(20, 10))\nax.imshow(viirs, vmin=0, vmax=30)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we see the total light emitted for each area:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"ntl_sums = zonal_stats(\n plants_buff,\n viirs,\n affine=viirs_aff,\n stats=[\"sum\"],\n nodata=-999,\n)\nplants[\"NTL\"] = [int(n[\"sum\"]) for n in ntl_sums]\nplants[[\"NTL\"]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As expected, high NTL from San Juan and Caguas, but much less from the other plant sites that aren't in major towns."},{"metadata":{},"cell_type":"markdown","source":"## Now let's check the wind roses for each site"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"for idx, row in plants.iterrows():\n print(\"Site:\", idx)\n x, y = row[\"xy\"]\n gfs_df = get_wind_data_point(gfs_path, x=x, y=y, disable_tqdm=True)\n rose = make_rose(gfs_df)\n fig, ax = draw_rose(rose, size=5)\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"They're all pretty similar, and all just about majority easterly wind, although the sites on the south coast do have a bit of a northerly component, and a bit weaker pverall. for simplicity let's just do all our wind profiles east-west. And let's set this decision into our table (remember 90° is east!)."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plants[\"wind_dir\"] = 90\nplants[[\"wind_dir\"]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Now the curve fit procedure for each site\nWhich gives the emissions (*E*) for each site in kg/hr."},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":false,"trusted":true},"cell_type":"code","source":"def get_E_value(x, y, sample_length=160):\n gfs_df = get_wind_data_point(gfs_path, x=x, y=y, disable_tqdm=True)\n \n s5p_arrs, aff = s5p_arr_timeseries(s5p_path, gfs_df, disable_tqdm=True)\n wind_cutoff = 2.5\n s5p_wind = [w[\"arr\"] for w in s5p_arrs if w[\"wind\"] >= wind_cutoff]\n s5p_wind_mean = np.nanmean(np.stack(s5p_wind), axis=0)\n \n line_df, sample_loc = get_line_density(s5p_wind_mean, s5p_aff, x, y, length=sample_length, width=40)\n xx = np.array(line_df.index.to_list())\n E, B, X, x0, std = fit_model(M, xx, line_df.val.to_list())\n yy = M(x=xx, E=E, B=B, X=X, x0=x0, std=std)\n plot_density(line_df, xx, yy, size=6)\n \n return E","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"plants[\"E_kghr\"] = 0\nfor idx, row in plants.iterrows():\n print(\"Site:\", idx)\n x, y = row[\"xy\"]\n E = get_E_value(x, y)\n plants.at[idx, \"E_kghr\"] = E * sec_per_hour * no2_mol_mass","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Clearly there is something wrong with the Costa Sur/EcoElectrica site! The line density profile makes no sense and there's no clear jump in NO2 behind the power plants. Either it is emitting very little NO2, or it just isn't showing up for some reason. Unfortunately I'll have to exclude it from our analysis for now.\n\nLet's have a look at the values for E in kg/hr for the remaining sites (remember this is the *total* E at each point, for electricity and other sources):"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plants = plants.drop([\"Costa Sur, EcoEléctrica\"])\nplants[[\"E_kghr\"]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Now we can solve for the unknown values\nWe have\n\n*E = E_elec + f_other x NTL*\n\nfor each of the points. And we know Caguas has *E_elec = 0* which let's us solve for *f_other*:\n\n*f_other = (E - E_elec)/NTL = E/NTL*"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"f_other = float(plants.at[\"Caguas\", \"E_kghr\"] / plants.at[\"Caguas\", \"NTL\"])\nprint(f\"f_other = {f_other:.4f}\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we can use *f_other* and *NTL* to calculate *Eelec* for each of the other stations"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plants[\"E_elec_kghr\"] = (plants[\"E_kghr\"] - f_other * plants[\"NTL\"]).astype(int)\nplants[[\"Capacity (MW)\", \"Hourly Generation (MWh)\", \"NTL\", \"E_kghr\", \"E_elec_kghr\"]]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And then we just use the estimated generation for each plant to calculate its emissions factor and then we can get a weighted emissions factor!"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plants[\"EF\"] = round(plants[\"E_elec_kghr\"] / plants[\"generation_hourly\"], 2)\nplants.at[\"Caguas\", \"EF\"] = 0\nplants[[\"Primary Fuel\", \"Capacity (MW)\", \"Hourly Generation (MWh)\", \"EF\"]]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"EF_mean = np.average(plants[\"EF\"], weights=plants[\"generation_hourly\"])\nprint(f\"Weighted average EF is {EF_mean:.2f} kg/MWh\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Is that it?? Are we done yet?\nJust about. We must bring back in the excluded power plants, some of which have emissions that we couldn't measure, and others which have no emissions."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"eia = pd.read_csv(eia_path)\neia[\"generation_hourly\"] = eia[\"Netgen_Total\"] / (265 * 24)\ntotal_gen = eia.generation_hourly.sum()\nfossil_codes = [\"NG\", \"BIT\", \"DFO\", \"RFO\"]\ntotal_fossil_gen = eia.loc[eia.Fuel_Code.isin(fossil_codes)].generation_hourly.sum()\ncaptured_gen = plants.generation_hourly.sum()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(f\"Total average hourly generation for Puerto Rico:\\t{total_gen:.0f} MWh\")\nprint(f\"Of which, the amount that is polluting:\\t\\t\\t{total_fossil_gen:.0f} MWh\\t({(100*total_fossil_gen/total_gen):.0f}% of total)\")\nprint(f\"Of which, the amount included in the analysis:\\t\\t{captured_gen:.0f} MWh\\t\\t({(100*captured_gen/total_gen):.0f}% of total)\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There's no single way to deal with this, but we can do our best. We can either:\n1. Assume our EF above applies to all energy produced on the island;\n2. Assume the rest have no emissions;\n3. Assume they follow general trends of pollution for their fuel type (coal, oil, or gas).\n\nNumber 1 is the number we already calculated, which is probably an overestimate, while number 2 will be significantly lower and is probably an under-estimate. Let's quickly try method 3 and see what it gives us.\n\nAll the emissions from the two San Juan plants are oil. Therefore we can take that site's EF as the Puerto Rico EF for oil. The other site we measured was a mix of coal and oil. Using the oil EF we calculate, we can use the weighted generation (by oil and by coal) to calculate the EF for coal. We have no information about gas... But luckily it doesn't account for much of Puerto Rico's generation. Going by our generic EF values from way up top, let's say it's five times lower than EF for oil. Thus we have for all three:"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"EF_oil = plants.at[\"San Juan CC, Palo Seco\", \"EF\"]\nprint(f\"EF for oil: {EF_oil:.2f} kg/MWh\")\n\ngen_coal = int(eia.loc[eia.Plant_Name == \"AES Puerto Rico\", \"generation_hourly\"])\ngen_oil = int(eia.groupby(\"Plant_Name\").sum().at[\"Aguirre Plant\", \"generation_hourly\"])\nEF_coal = (((plants.at[\"A.E.S. Corp., Aguirre\", \"EF\"]) * (gen_coal + gen_oil)) - (gen_oil * EF_oil))/gen_coal\nprint(f\"EF for coal: {EF_coal:.2f} kg/MWh\")\n\nEF_gas = EF_oil / 5\nprint(f\"EF for gas: {EF_gas:.2f} kg/MWh\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we're ready to extrapolate our total emissions."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def calc_emissions(row):\n if row.Fuel_Code in [\"NG\"]:\n return EF_gas * row[\"generation_hourly\"]\n elif row.Fuel_Code in [\"BIT\"]:\n return EF_coal * row[\"generation_hourly\"]\n elif row.Fuel_Code in [\"DFO\", \"RFO\"]:\n return EF_oil * row[\"generation_hourly\"]\n else:\n return 0\neia[\"emissions_kghr\"] = eia.apply(calc_emissions, axis=\"columns\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"EF_method3 = eia.emissions_kghr.sum() / eia.generation_hourly.sum()\nprint(f\"EF if we extrapolate to other plants: {EF_method3:.2f} kg/MWh\")\nprint(\"This is probably the best estimate!\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## So there we have it, final answer!\nOr a bit lower if we assume the other plants don't pollute, or a bit higher if we assume uniform emissions.\n\n# Bonus Points\n## We can now calculate the marginal emissions factor\nFor a marginal emissions factor to be meaningful, we'd need a lot more generation data, down to daily or better resolution. This would enable an understanding of, at each point in time, what the environmental impact of a decision to use electricity would be. Even an average factor would need to take into account which plants are historically at the margin of additional power demand. Since we don't have that, let's just assume that marginal power comes from fossil fuel plants. There are quite a few ways we could go about this in the absence of data...\n\nLet's say that coal is unlikely to be on the margin (oil and gas are hapipier being ramped up and down), so our marginal power will come from either gas or oil. Then we can calculate aan average (weighted by generation) of the emission factors of each of these fuels."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"gen_gas = eia.loc[eia.Fuel_Code == \"NG\"].Netgen_Total.sum()\ngen_oil = eia.loc[eia.Fuel_Code.isin([\"RFO\", \"DFO\"])].Netgen_Total.sum()\nEF_marginal = (EF_gas * gen_gas + EF_oil * gen_oil)/(gen_gas + gen_oil)\nprint(f\"Marginal EF: {EF_marginal:.2f} kg/MWh\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## A bit lower, but just because we left out coal\nThis really requires a bit more knowledge about Puerto Rico's daily generation and price-demand curves!\n\n## What about calculating the EF on a monthly basis?\nIn theory easy (just apply the method to only a month of data at a time), but as there's only one NO2 image per day, we start to stretch it a bit thin. Will do it and we can see what the results look like and if they correlate with anything. No expectation that they will change dramatically, as renewable (wind, solar, hydro) penetration isn't high enough (~6%) to drastically change the proportion of electricity coming from fossil (polluting) sources.\n\nLet's see how we would do this just for San Juan. First we create a new dict with separate NO2 and wind data for each month."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def calc_monthly(x, y):\n gfs_df = get_wind_data_point(gfs_path, x=x, y=y, disable_tqdm=True)\n s5p_arrs, aff = s5p_arr_timeseries(s5p_path, gfs_df, disable_tqdm=True)\n \n s5p_arrs = sorted(s5p_arrs, key=lambda x: x[\"date\"])\n s5p_monthly = {}\n for d in s5p_arrs:\n year = d[\"date\"].year\n month = d[\"date\"].month\n ym = f\"{year}-{month:02}\"\n if ym not in s5p_monthly.keys():\n s5p_monthly[ym] = []\n s5p_monthly[ym].append(d)\n \n results = {}\n for month, vals in s5p_monthly.items():\n wind_cutoff = 2.5\n s5p_wind = [w[\"arr\"] for w in vals if w[\"wind\"] >= wind_cutoff]\n s5p_wind_mean = np.nanmean(np.stack(s5p_wind), axis=0)\n line_df, sample_loc = get_line_density(s5p_wind_mean, aff, x, y, length=160, width=40)\n xx = np.array(line_df.index.to_list())\n try:\n E, B, X, x0, std = fit_model(M, xx, line_df.val.to_list())\n except RuntimeError:\n E = np.nan\n results[month] = int(E * sec_per_hour * no2_mol_mass)\n return results","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"x, y = plants.at[\"San Juan CC, Palo Seco\", \"xy\"]\nE_monthly = calc_monthly(x, y)\nprint(\"Monthly emissions for San Juan area (from all sources!)\")\nmonthly = pd.DataFrame(E_monthly.values(), index=E_monthly.keys(), columns=[\"E_kghr\"])\nmonthly.index = pd.to_datetime(monthly.index)\nmonthly.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"At this point we could repeat the above procedure, using NTL each month (NTL data is available as monthly composites!) to disaggregate power and non-power emissions, and using the monthly generation from each plant to determine a monthly EF factor. However, I think this is sufficient as a demonstration of the potential of this method, Anyway I think we might start stretching our assumptions a bit thin to keep going without stopping for a validity check on how we've done so far!\n\nTo get a quick answer out of this, let's just divide the emissions each month by Puerto Rico's total generation for that month (remembering that the generation is from 2018, whereas our NO2 data is from 2018-2019, but we assume each month has similar characteristics regardless of year). Then we'll just normalise the results to have a mean equal to the annual EF we already calculated!"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"monthly[\"EF_high\"] = np.nan\nfor idx, row in monthly.iterrows():\n month = idx.strftime(\"%B\")\n monthly.at[idx, \"EF_high\"] = round(row[\"E_kghr\"] / (eia[f\"Netgen_{month}\"].sum() / (24 * 365)), 2)\nmonthly[\"EF\"] = (monthly[\"EF_high\"] * EF_method3 / monthly[\"EF_high\"].mean()).round(2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And this is what it looks like!"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def plot_monthly_EF(df):\n sns.set_style(\"whitegrid\")\n fig, ax = plt.subplots(figsize=(16,6))\n df[[\"EF\"]].plot(ax=ax, linewidth=4)\n ax.set_ylim([0, 3])\n ax.set_title(\"Emissions Factor (EF) by date\", fontsize=16)\n ax.set_ylabel(\"EF (kg/MWh)\", fontsize=16)\n ax.set_xlabel(\"Date\", fontsize=16)\n ax.tick_params(axis='both', which='both', labelsize=12)\n \n ax.xaxis.grid(False)\n ax.legend()\n sns.despine(left=True)\n plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plot_monthly_EF(monthly)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## And now we're seriously done!\n\n# Caveats and Recommendations\nAs mentioned, I decided to keep it as simple as possible for this demonstration. There is also quite a long chain of assumptions that could build up to several inaccuracies.\n\nThe most obvious improvement would be to include a larger patch of wind and NO2 data, extending much further to each side of Puerto Rico (especially downwind!). This should give much more freedom in creating our wind profiles. As mentioned, we should also use atmospheric data that matches month for month with our generation data, but this isn't always possible if we want the latest numbers! (I haven't managed to find 2019 generation numbers for Puerto Rico, but they probably exist at the end of the right email address).\n\nOne big assumption was assuming non-power pollution scales linearly with night time lights (NTL). We should also investigate other data sources such as population distribution, road networks etc, and explore the possibility of a nonlinear relationship (e.g. more night time lights in one area indicates more intensive industry, which could create a much more than linear increase in NO2 emissions). Our choice of NTL sampling area (circle with 5 km radius) was also essentially random and needs to be investigated.\n\nWe were lucky that Puerto Rico's wind is so well behaved. A more general analysis would require rotating the NO2 fields (as has been done in the literature) or separatingly sampling along major wind directions (which I wrote the code for above but didn't use). The line density profile also assumes a near point source of emissions, but we know that isn't the case. We might need to investigate including multiple model functions to fit one density line, each representing a different emission source (this might allow us to model the plants further west under San Juan's plume).\n\nWe should also look into additional effects on NO2, such as temperature, humidity, precipitation etc (in this analysis I haven't even touched the [GLDAS](https://developers.google.com/earth-engine/datasets/catalog/NASA_GLDAS_V021_NOAH_G025_T3H) dataset). However, without decent training data it will be difficult to include these other data sources that aren't as easy to analyse with a simple, logical approach.\n\nTo apply this method in other areas, there are a few things which would need to be considered (in addition to those mentioned above). The assumption about no wind blowing in would generally not be true, though this shouldn't affect the line density procedure much unless the plumes are very close. For areas with more variable wind, the analysis would either have to be done separately for different wind directions, or the NO2 fields rotated according to each day's wind. The process of disaggregating power and non-power pollution would need to be assessed, and in dense areas it is likely that many plants would be under the plumes of others, in which case the suggestion above of multiple model functions should be explored."},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}