{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Visualising Wind Effects\n\nWhen the wind is a blowin' the emissions be a showin', or at least, that's what I hope. NO2 emissions are a gas, which means they are easily influenced by wind. Prior research shows that wind can significantly effect the dispersion profile of NO2 gas so we should take this into account when trying to estimate the emissions from the power generation facilities."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import os\n\nimport rasterio as rio\nfrom rasterio.plot import show, show_hist\nimport matplotlib.pyplot as plt\nfrom matplotlib import animation\nimport seaborn as sns\n\nfrom datetime import datetime\n\nimport numpy as np\nimport pandas as pd","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Read Data\n\nGet all data for NO2 and weather."},{"metadata":{"code_folding":[0,31],"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def load_sp5():\n    # Get the data filenames\n    no2_path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/s5p_no2/'\n    no2_files = [no2_path + f for f in os.listdir(no2_path)]\n    \n    data = []\n    \n    print(f'Reading {len(no2_files)} files')\n    for f in no2_files:\n        raster = rio.open(f)\n        data += [{\n            'tif': raster,\n            'filename': f.split('/')[-1],\n            'path': no2_path + f.split('/')[-1],\n            'measurement': 'no2',\n            'band1_mean': np.nanmean(raster.read(1)),\n            **raster.meta\n        }]\n        raster.close()\n        \n    # Get dates\n    for d in data:\n        d.update({'datetime': datetime.strptime(d['filename'][:23], 's5p_no2_%Y%m%dT%H%M%S')})\n\n    for d in data:\n        d['date'] = d['datetime'].date()\n        d['hour'] = int(datetime.strftime(d['datetime'], '%H'))\n        d['weekday'] = d['datetime'].weekday()  # Mon = 0\n\n    return data\n\ndef load_gfs():\n    # Get the data filenames\n    no2_path = '/kaggle/input/ds4g-environmental-insights-explorer/eie_data/gfs/'\n    no2_files = [no2_path + f for f in os.listdir(no2_path)]\n    \n    data = []\n    \n    print(f'Reading {len(no2_files)} files')\n    for f in no2_files:\n        raster = rio.open(f)\n        data += [{\n            'tif': raster,\n            'filename': f.split('/')[-1],\n            'path': no2_path + f.split('/')[-1],\n            'measurement': 'weather',\n            'band1_mean': np.nanmean(raster.read(1)),\n            **raster.meta\n        }]\n        raster.close()\n        \n    # Get dates\n    for d in data:\n        d.update({'datetime': datetime.strptime(d['filename'], 'gfs_%Y%m%d%H.tif')})\n\n    for d in data:\n        d['date'] = d['datetime'].date()\n        d['hour'] = int(datetime.strftime(d['datetime'], '%H'))\n        d['weekday'] = d['datetime'].weekday()  # Mon = 0\n\n    return data\n\ndata = load_sp5() + load_gfs()\n# data.append(load_gfs)\n","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"data[0]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(data)\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Weather\n\nWe'll discuss wind, but first a mention on reflectance.\n\n### High Reflectance\n\nWhen there is high reflectance, such as cloud cover or snow, the satellite can have a hard time measuring a good signal, as such these values should be ignored. The dataset in use here has already been processed to remove cases of significant cloud coverage but we can still take note of the cloud values as they may affect our readings.\n\nThankfully, Puerto Rico is a warm state, meaning that it doesn't receive any snow throughout the year.\n\nReading:\n\n- https://traveltips.usatoday.com/weather-climate-puerto-rico-50065.html\n\n### Wind\n\nStrong winds move the air in the atmosphere which means that the spatial accuracy of measurements quickly becomes incorrect. Thankfully, these measurements can be corrected for if we have a good understanding of the wind. In general, it is sensible to look at days with less wind.\n\nFor a more accurate analysis, we can remove from our dataset dates of measurements which correspond to strong wind conditions. This will more accurately allocate the NO2 readings to the spatial source.\n\nThe wind measurements occur every 6 hours, but the NO2 data is once per day. Since we're trying to move any large cases of wind, we can aggregate the wind values to daily, using the mean/max. This can then be joined with the NO2 data.\n\nOther things to consider:\n\n- Counter winds could be OK, e.g. strong easterly followed by strong westerly wind.\n- Wind fluctuates a lot in time, so 6 hourly measurements may not be enough to detect the impact on emissions.\n\n\nLet's calculate the magnitude's and have a look."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Add total wind magnitude\nfor d in data:\n    if d['measurement'] == 'weather':\n        r = rio.open(d['path'])\n        d.update({\n            'avg_wind_magnitude': \n            np.nanmean(np.sqrt(np.square(r.read(4)) + np.square(r.read(5)))),\n            'max_wind_magnitude': \n            np.nanmax(np.sqrt(np.square(r.read(4)) + np.square(r.read(5))))\n        })\n        r.close()\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.hist([d['max_wind_magnitude'] for d in data if d['measurement'] == 'weather'])\nplt.hist([d['avg_wind_magnitude'] for d in data if d['measurement'] == 'weather'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see there are some readings with a much higher mean (6+), but otherwise the data seems to be distributed with some positive skew.\n\nThe units here are m/s. The average magnitude is 4m/s or 14km/h. For reference on the scale, Puerto Rico is approximately 50km from top to bottom and 150km from left to right. So we see that these winds can definitely move the emissions around on this scale.\n\nTo continue, we can filter out the lower quartile and compare it with the values in the upper quartile.\n\nLet's plot the daily values to see any patterns, with and without the high wind days."},{"metadata":{"code_folding":[0],"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def raster_median(rasters, band=1):\n    tmp = []\n    for r in rasters:\n        r = rio.open(r.name)\n        tmp += [r.read(band)]\n        r.close()\n    return [np.nanmedian(np.stack(tmp), axis=(0))]","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(data)\ndf_daily = df.groupby('date').agg({\n    'avg_wind_magnitude': 'median',\n    'max_wind_magnitude': 'max',\n    'tif': raster_median,\n    'measurement': 'first'\n}).rename(columns={'tif': 'values'})\ndf_daily.plot(kind='line')","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"quartile25 = df_daily['avg_wind_magnitude'].quantile(0.25)\nquartile75 = df_daily['avg_wind_magnitude'].quantile(0.75)\n\n# df_daily.plot(kind='line')\nprint('Low Wind Days')\ndf_daily.query('avg_wind_magnitude < @quartile25').plot(kind='line', marker='.', linestyle='')\nprint(\"High Wind Days\")\ndf_daily.query('avg_wind_magnitude > @quartile75').plot(kind='line', marker='.', linestyle='')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Visualising\n\nWe can now compare the difference in the emissions data from removing the wind. How does our data compare?\n\nWe'll use the methods described here: https://www.kaggle.com/jyesawtellrickson/averaging-satellite-data."},{"metadata":{"code_folding":[0],"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def plot_average_raster(rasters, band=1, output_file='tmp.tif', avg='mean'):\n    all_no2s = []\n    print(f'Processing {len(rasters)} files')\n    for r in rasters:\n        if r.closed:\n            r = rio.open(r.name)\n        all_no2s += [r.read()[band-1, :, :]]\n        r.close()\n    temporal_no2 = np.stack(all_no2s)\n    print(\"Done\")\n    \n    if avg == 'mean':\n        avg_no2 = np.nanmean(temporal_no2, axis=(0))\n    else:\n        avg_no2 = np.nanmedian(temporal_no2, axis=(0))\n\n    raster = rasters[0]\n    \n    new_dataset = rio.open(\n        output_file,\n        'w',\n        driver=raster.driver,\n        height=raster.height,\n        width=raster.width,\n        count=1,\n        dtype=avg_no2.dtype,\n        crs=raster.crs,\n        transform=raster.transform,\n    )\n    \n    new_dataset.write(avg_no2*10**5, 1)\n    new_dataset.close()\n    \n    tmp = rio.open(output_file)\n    \n    min_val = np.nanmin(tmp.read(1))\n    max_val = np.nanmax(tmp.read(1))\n    print('Ranges from {:.2E} to {:.2E}'.format(min_val, max_val))\n    \n    # https://rasterio.readthedocs.io/en/latest/topics/plotting.html\n    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 8))\n\n    # Augment the data so that it can be plotted nicely\n    # mult_fact = 10**(-round(np.log10(min_val))+1)\n    #data = tmp.read(1)*10**mult_fact\n    \n    show(tmp, transform=tmp.transform, ax=ax1)\n    \n    show(tmp, cmap='Greys_r', interpolation='none', ax=ax2)\n    show(tmp, contour=True, ax=ax2)\n\n    plt.show()\n    \n    return tmp","execution_count":null,"outputs":[]},{"metadata":{"code_folding":[0],"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def plot_diff_raster(rasters_1, rasters_2, band=1, output_file='tmp.tif'):\n    all_no2s_1 = []\n    all_no2s_2 = []\n    # Read raster 1\n    print(f'Processing {len(rasters_1)} files')\n    for r in rasters_1:\n        r = rio.open(r.name)\n        all_no2s_1 += [r.read()[band-1, :, :]]\n        r.close()\n    temporal_no2_1 = np.stack(all_no2s_1)\n    # Read raster 2\n    print(f'Processing {len(rasters_2)} files')\n    for r in rasters_2:\n        r = rio.open(r.name)\n        all_no2s_2 += [r.read()[band-1, :, :]]\n        r.close()\n    temporal_no2_2 = np.stack(all_no2s_2)\n        \n    # Calculate averages\n    avg_no2_1 = np.nanmean(temporal_no2_1, axis=(0))\n    avg_no2_2 = np.nanmean(temporal_no2_2, axis=(0))\n\n    avg_no2 = avg_no2_2 - avg_no2_1\n    \n    raster = rasters_1[0]\n    \n    new_dataset = rio.open(\n        output_file,\n        'w',\n        driver=raster.driver,\n        height=raster.height,\n        width=raster.width,\n        count=1,\n        dtype=avg_no2.dtype,\n        crs=raster.crs,\n        transform=raster.transform,\n    )\n    \n    new_dataset.write(avg_no2*10**5, 1)\n    new_dataset.close()\n    \n    tmp = rio.open(output_file)\n    \n    print('Ranges from {:.2E} to {:.2E}'.format(np.nanmin(tmp.read(1)),np.nanmax(tmp.read(1))))\n    \n    # https://rasterio.readthedocs.io/en/latest/topics/plotting.html\n    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 8))\n    \n    show(tmp, transform=tmp.transform, ax=ax1)\n    \n    show(tmp, cmap='Greys_r', interpolation='none', ax=ax2)\n    show(tmp, contour=True, ax=ax2)\n\n    plt.show()\n    \n    return tmp","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\"Low Wind Days\")\nplot_average_raster(\n    df[df.date.isin(\n        df_daily.query('avg_wind_magnitude < @quartile25*0.9'\n                      ).index.tolist())].query('measurement == \"no2\"').tif.tolist()\n)\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\"High Wind Days\")\navg = plot_average_raster(\n    df[df.date.isin(\n        df_daily.query('avg_wind_magnitude > @quartile75'\n                      ).index.tolist())].query('measurement == \"no2\"').tif.tolist()\n)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It's also constructive to plot the difference between the days with less wind vs. the days with a lot of wind."},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_diff_raster(\n    df[df.date.isin(\n        df_daily.query('avg_wind_magnitude < @quartile25*0.9'\n                      ).index.tolist())].query('measurement == \"no2\"').tif.tolist(),\n    df[df.date.isin(\n        df_daily.query('avg_wind_magnitude > @quartile75'\n                      ).index.tolist())].query('measurement == \"no2\"').tif.tolist()    \n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So we can see quite well that on the windy days, there is less contrast, and the clump of pollution around San Juan is more spread out to the west.\n\nWe also see this in the values on the contour plot. On the days without wind, the concentration is stronger (inner circle of 6.8) compared to the days with wind (inner circle of 6.0).\n\nThe difference plot shows a fairly even spread outside of the big decrease in San Juan.\n\nThere's a potential that the date / time of our reading could be leading to the difference in levels, so let's check if there's a correlation between any of the variables with our wind."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"corr_no2 = df.query('measurement == \"no2\"').corr()\ncorr_qgs = df.query('measurement == \"weather\"').corr()\n\n\n# https://rasterio.readthedocs.io/en/latest/topics/plotting.html\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n\nsns.heatmap(\n    corr_no2, \n    vmin=-1, vmax=1, center=0,\n    cmap=sns.diverging_palette(20, 220, n=200),\n    square=True,\n    ax=ax1\n)\nax1.set_xticklabels(\n    ax1.get_xticklabels(),\n    rotation=45,\n    horizontalalignment='right'\n);\nplt.title(\"NO2\")\n\nsns.heatmap(\n    corr_qgs, \n    vmin=-1, vmax=1, center=0,\n    cmap=sns.diverging_palette(20, 220, n=200),\n    square=True,\n    ax=ax2\n)\nax2.set_xticklabels(\n    ax2.get_xticklabels(),\n    rotation=45,\n    horizontalalignment='right'\n);\nplt.title('Weather')\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We do see some correlations here.\n\nFor the left plot, the band 1 mean (NO2 concentration) is weakly correlated with the hour of the day. This makes sense, since the emissions are largely driven by cooling in the afternoon and appliance usage in the evening. We should check that the avg. hour for the two plots is roughly the same.\n\nWe can also see that the hour is weakly correlated with the avg wind magnitude, however this is not considered to be a problem as the values for the wind measurements are 6-hourly and averaged out for the day."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print(\"Less wind values, avg. hour: {}\".format(\ndf[df.date.isin(\n    df_daily.query('avg_wind_magnitude < @quartile25*0.9'\n                  ).index.tolist())].query('measurement == \"no2\"').hour.mean()\n))\n\nprint(\"More wind values, avg. hour: {}\".format(\ndf[df.date.isin(\n    df_daily.query('avg_wind_magnitude > @quartile75'\n                  ).index.tolist())].query('measurement == \"no2\"').hour.mean()\n))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The values are very close, so it seems we're OK.\n\n## Conclusions\n\nThe wind can have a significant effect on the emissions distributions as shown by removing the readings with large amounts of wind.\n\nIn order to get more precise emissions readings it's better to use days without wind.\n\nNext steps:\n\n- There's research which shows the windy days can be utilised, through some smart calculations, this could be applied here in order to make more use of the readings and apply the analytical method to locations that have large amounts of wind (more universal methodology).\n- How long does NO2 last in the readings? If there were strong winds the hour / 6 hours / day before, would the NO2 from a source be lingering in the region? Should we exclude a certain window around strong wind events?\n\n\nReading:\n- <https://www.atmos-chem-phys.net/16/5283/2016/acp-16-5283-2016.pdf>"}],"metadata":{"hide_input":false,"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.9"},"toc":{"base_numbering":1,"nav_menu":{},"number_sections":true,"sideBar":true,"skip_h1_title":true,"title_cell":"Table of Contents","title_sidebar":"Contents","toc_cell":false,"toc_position":{},"toc_section_display":true,"toc_window_display":true},"varInspector":{"cols":{"lenName":16,"lenType":16,"lenVar":40},"kernels_config":{"python":{"delete_cmd_postfix":"","delete_cmd_prefix":"del ","library":"var_list.py","varRefreshCmd":"print(var_dic_list())"},"r":{"delete_cmd_postfix":") ","delete_cmd_prefix":"rm(","library":"var_list.r","varRefreshCmd":"cat(var_dic_list()) "}},"types_to_exclude":["module","function","builtin_function_or_method","instance","_Feature"],"window_display":false}},"nbformat":4,"nbformat_minor":1}