{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Averaging Satellite Data\n\nSatellite data doesn't have simple, consistent measures of locations to work with. Due to the variability in the location of the satellite and various other conditions, the exact locations vary. This makes getting a measurement for a specific point more difficult.\n\nIn this notebook we'll explore the Sentinel-5P dataset measurement locations to see if we can make an average."},{"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\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\nWe'll work with the Sentinel-5P dataset.\n\nPrior work: https://www.kaggle.com/jyesawtellrickson/understanding-sentinel-5p-offl-no2"},{"metadata":{"_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            'measurement': 'no2',\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'] = datetime.strftime(d['datetime'], '%H')\n        d['weekday'] = d['datetime'].weekday()  # Mon = 0\n\n    return data\n\ndata = load_sp5()","execution_count":null,"outputs":[]},{"metadata":{"_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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Affine Transformations\n\nThe measurements from satellite are recorded with two important location information. These are:\n\n- CRS: The coordinate reference system that the measurements are referenced to. This is the world situation.\n- Transform: The affine transformation that used to get from the row/column indexes stored to x/y values in the coordinate reference system.\n\nWith the above two pieces of information we can take the values stored in the TIFF file and transform them to values in an x/y plane in our given coordinate reference system. In this case the CRS is EPSG:4326 for all the readings.\n\nFurther reading:\n- Affine Transformations: https://stackabuse.com/affine-image-transformations-in-python-with-numpy-pillow-and-opencv/\n- EPSG:4326: https://epsg.io/4326\n\nLet's start by investigating the affine transformations being used."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# Get the affine transformation values\naff = pd.DataFrame(df['transform'].values.tolist())\naff.head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"for col in aff.columns:\n    print(aff[col].value_counts(), '\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The letters a through i refer to the values in the affine transformation matrix:\n\n- a/e are responsible for scaling\n- b/d are responsible for shearing\n- c/f are responsible for translations\n- g/h/i are set to 0, 0 and 1 respectively\n\nHere we see that the scaling in both directions is the same, 0.004492. This number also represents the width of one of our readings, the pixel width. We have each reading equal to 0.004492 degrees, roughly 500m. The shearing is 0. The translations are varied, but are roughly -67.325 in the x and 17.901 in the y. The variance seen is no more than 0.002, which is less than a half of the pixel width.\n\nSo we see that in most cases are coordinate systems are aligned, but not all. When they do vary, it's quite a small variance, which means as a first attempt, taking a straight average should give a decent result. To be more accurate, we could ignore the readings that show 0.002 variance in the x direction as they may cause some more sizable overlap and their quantity is small. \n\nReading:\n- https://stackabuse.com/affine-image-transformations-in-python-with-numpy-pillow-and-opencv/\n- https://en.wikipedia.org/wiki/Decimal_degrees\n- http://www.csgnetwork.com/degreelenllavcalc.html"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"df = pd.merge(df, aff, left_index=True, right_index=True, how='inner')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Creating the Average\n\nLet's try taking the average and plotting. We can compare using the values that have varying translations or the ones without. We can also compare taking mean vs. median."},{"metadata":{"code_folding":[],"_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    \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, 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)),\n                                                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=(6, 6))\n\n    show(tmp, transform=tmp.transform, ax=ax1)\n    \n    show((tmp, 1), cmap='Greys_r', interpolation='none', ax=ax2)\n    show((tmp, 1), 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('All files, mean', '\\n')\n\ntmp = plot_average_raster(df.tif.tolist(), avg='mean')","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('No translation files, mean', '\\n')\n\ntmp = plot_average_raster(df.query('f == 17.901140352016327 and c == -67.32431391288841').tif.tolist(),\n                          avg='mean')","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"print('All files', '\\n')\n\ntmp = plot_average_raster(df.tif.tolist(), avg='median')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There's very little difference made by removing the translated files, so it should be fine to proceed as is. Similarly, mean and median show little differences in the results, though the boundaries in the contour plot seem to be less smooth.\n\n\nAlso, for some reason our plots are upside down, but we can definitely see strong effects here. Namely, the capital of Puerto Rico, San Juan shows the peak emissions with those emissions also going out to the left strongly (maybe this is wind, or another city). We also see some other isolated areas of high / low concentrations towards the centre-east of the island.\n\n## Conclusions\n\nWhile averages are great because they allow us to get a quick picture of what is going on, they smooth out all the interesting effects that be required to capture a better idea about the emissions effects of power plants.\n\nNext up we should look deeper into:\n - wind strength\n - day of the week\n - time of day\n - special events"}],"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}