Slope_Stats_FINAL.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Slope_Stats_FINAL.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyNl7qpaTRe8msDrCoujBtlk",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Ario20/Melite/blob/main/Slope_Stats_FINAL.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Version: 1.1\n",
        "Author: Kenneth Cassar\n",
        "Date: 10-07-2022\n",
        "\n",
        "The following Python script was adapted from Vollrath, Mullissa and Reiche 2020: \"Angular-based radiometric slope correction for Sentinel-1 on Google Earth Engine\". The GitHub repository for the paper is: https://github.com/ESA-PhiLab/radiometric-slope-correction\n",
        "\n",
        "This script implements the slope correction on the roi and calculates various statistics used in the evaluation of the best slope correction model according to the land cover type.\n",
        "\n",
        "The script is implemented in Google Colab and statistiscal data are saved in the user's Google Drive"
      ],
      "metadata": {
        "id": "128NxQrNohaq"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eAjuYfA3q0s7",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "7de69682-1a9f-4c6b-c2c9-92088029f07b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions. If the web browser does not start automatically, please manually browse the URL below.\n",
            "\n",
            "    https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=YMfU9QK7Mxqd0yZdJoKSejat_-EMg24KvDrYxm7Ob1w&tc=vN4NL1mws-CsBaYNItiFwbR5yXJulOcDoak0u4PXHeM&cc=53aCjg3bdL6PZtAvcY_NwwYDftxxsaf0X15r88unQXg\n",
            "\n",
            "The authorization workflow will generate a code, which you should paste in the box below. \n",
            "Enter verification code: 4/1AdQt8qgfYo56MRHk_IglfWXXSOTVtnUsELj0BGji1FGqkPibVEmE5fO4yS0\n",
            "\n",
            "Successfully saved authorization token.\n"
          ]
        }
      ],
      "source": [
        "import ee\n",
        "# Trigger the authentication flow.\n",
        "ee.Authenticate()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "ee.Initialize()"
      ],
      "metadata": {
        "id": "p_6f7TOGrBpL"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "1. Defining the slope correction module"
      ],
      "metadata": {
        "id": "1zpy-hSVqfyg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def slope_correction(collection, elevation, model, buffer=0):\n",
        "    '''This function applies the slope correction on a collection of Sentinel-1 data\n",
        "       \n",
        "       :param collection: ee.Collection of Sentinel-1\n",
        "       :param elevation: ee.Image of DEM\n",
        "       :param model: model to be applied (volume/surface)\n",
        "       :param buffer: buffer in meters for layover/shadow amsk\n",
        "        \n",
        "        :returns: ee.Image\n",
        "    '''\n",
        "    \n",
        "    def _volumetric_model_SCF(theta_iRad, alpha_rRad):\n",
        "        '''Code for calculation of volumetric model SCF\n",
        "        \n",
        "        :param theta_iRad: ee.Image of incidence angle in radians\n",
        "        :param alpha_rRad: ee.Image of slope steepness in range\n",
        "        \n",
        "        :returns: ee.Image\n",
        "        '''\n",
        "        \n",
        "        # create a 90 degree image in radians\n",
        "        ninetyRad = ee.Image.constant(90).multiply(np.pi/180)\n",
        "        \n",
        "        # model\n",
        "        nominator = (ninetyRad.subtract(theta_iRad).add(alpha_rRad)).tan()\n",
        "        denominator = (ninetyRad.subtract(theta_iRad)).tan()\n",
        "        return nominator.divide(denominator) \n",
        "    \n",
        "    \n",
        "    def _surface_model_SCF(theta_iRad, alpha_rRad, alpha_azRad):\n",
        "        '''Code for calculation of direct model SCF\n",
        "        \n",
        "        :param theta_iRad: ee.Image of incidence angle in radians\n",
        "        :param alpha_rRad: ee.Image of slope steepness in range\n",
        "        :param alpha_azRad: ee.Image of slope steepness in azimuth\n",
        "        \n",
        "        :returns: ee.Image\n",
        "        '''\n",
        "        \n",
        "        # create a 90 degree image in radians\n",
        "        ninetyRad = ee.Image.constant(90).multiply(np.pi/180)\n",
        "        \n",
        "        # model  \n",
        "        nominator = (ninetyRad.subtract(theta_iRad)).cos()\n",
        "        denominator = (alpha_azRad.cos()\n",
        "          .multiply((ninetyRad.subtract(theta_iRad).add(alpha_rRad)).cos()))\n",
        "\n",
        "        return nominator.divide(denominator)\n",
        "\n",
        "\n",
        "    def _erode(image, distance):\n",
        "      '''Buffer function for raster\n",
        "\n",
        "      :param image: ee.Image that shoudl be buffered\n",
        "      :param distance: distance of buffer in meters\n",
        "        \n",
        "      :returns: ee.Image\n",
        "      '''\n",
        "      \n",
        "      d = (image.Not().unmask(1)\n",
        "          .fastDistanceTransform(30).sqrt()\n",
        "          .multiply(ee.Image.pixelArea().sqrt()))\n",
        "    \n",
        "      return image.updateMask(d.gt(distance))\n",
        "    \n",
        "    \n",
        "    def _masking(alpha_rRad, theta_iRad, buffer):\n",
        "        '''Masking of layover and shadow\n",
        "        \n",
        "        \n",
        "        :param alpha_rRad: ee.Image of slope steepness in range\n",
        "        :param theta_iRad: ee.Image of incidence angle in radians\n",
        "        :param buffer: buffer in meters\n",
        "        \n",
        "        :returns: ee.Image\n",
        "        '''\n",
        "        # layover, where slope > radar viewing angle \n",
        "        layover = alpha_rRad.lt(theta_iRad).rename('layover')\n",
        "\n",
        "        # shadow \n",
        "        ninetyRad = ee.Image.constant(90).multiply(np.pi/180)\n",
        "        shadow = alpha_rRad.gt(ee.Image.constant(-1).multiply(ninetyRad.subtract(theta_iRad))).rename('shadow')\n",
        "        \n",
        "        # add buffer to layover and shadow\n",
        "        if buffer > 0:\n",
        "            layover = _erode(layover, buffer)   \n",
        "            shadow = _erode(shadow, buffer)  \n",
        "\n",
        "        # combine layover and shadow\n",
        "        no_data_mask = layover.And(shadow).rename('no_data_mask')\n",
        "        \n",
        "        return layover.addBands(shadow).addBands(no_data_mask)\n",
        "                        \n",
        "        \n",
        "    def _correct(image):\n",
        "        '''This function applies the slope correction and adds layover and shadow masks\n",
        "        \n",
        "        '''\n",
        "        \n",
        "        # get the image geometry and projection\n",
        "        geom = image.geometry()\n",
        "        proj = image.select(1).projection()\n",
        "        \n",
        "        # calculate the look direction\n",
        "        heading = (ee.Terrain.aspect(image.select('angle'))\n",
        "                                     .reduceRegion(ee.Reducer.mean(), geom, 1000)\n",
        "                                     .get('aspect'))\n",
        "                   \n",
        "\n",
        "        # Sigma0 to Power of input image\n",
        "        sigma0Pow = ee.Image.constant(10).pow(image.divide(10.0))\n",
        "\n",
        "        # the numbering follows the article chapters\n",
        "        # 2.1.1 Radar geometry \n",
        "        theta_iRad = image.select('angle').multiply(np.pi/180)\n",
        "        phi_iRad = ee.Image.constant(heading).multiply(np.pi/180)\n",
        "        \n",
        "        # 2.1.2 Terrain geometry\n",
        "        alpha_sRad = ee.Terrain.slope(elevation).select('slope').multiply(np.pi/180).setDefaultProjection(proj).clip(geom)\n",
        "        phi_sRad = ee.Terrain.aspect(elevation).select('aspect').multiply(np.pi/180).setDefaultProjection(proj).clip(geom)\n",
        "        \n",
        "        # we get the height, for export \n",
        "        height = elevation.setDefaultProjection(proj).clip(geom)\n",
        "        \n",
        "        # 2.1.3 Model geometry\n",
        "        #reduce to 3 angle\n",
        "        phi_rRad = phi_iRad.subtract(phi_sRad)\n",
        "\n",
        "        # slope steepness in range (eq. 2)\n",
        "        alpha_rRad = (alpha_sRad.tan().multiply(phi_rRad.cos())).atan()\n",
        "\n",
        "        # slope steepness in azimuth (eq 3)\n",
        "        alpha_azRad = (alpha_sRad.tan().multiply(phi_rRad.sin())).atan()\n",
        "\n",
        "        # local incidence angle (eq. 4)\n",
        "        theta_liaRad = (alpha_azRad.cos().multiply((theta_iRad.subtract(alpha_rRad)).cos())).acos()\n",
        "        theta_liaDeg = theta_liaRad.multiply(180/np.pi)\n",
        "\n",
        "        # 2.2 \n",
        "        # Gamma_nought\n",
        "        gamma0 = sigma0Pow.divide(theta_iRad.cos())\n",
        "        gamma0dB = ee.Image.constant(10).multiply(gamma0.log10()).select(['VV', 'VH'], ['VV_gamma0', 'VH_gamma0'])\n",
        "        ratio_gamma = (gamma0dB.select('VV_gamma0')\n",
        "                        .subtract(gamma0dB.select('VH_gamma0'))\n",
        "                        .rename('ratio_gamma0'))\n",
        "\n",
        "        if model == 'volume':\n",
        "            scf = _volumetric_model_SCF(theta_iRad, alpha_rRad)\n",
        "\n",
        "        if model == 'surface':\n",
        "            scf = _surface_model_SCF(theta_iRad, alpha_rRad, alpha_azRad)\n",
        "\n",
        "        # apply model for Gamm0_f\n",
        "        gamma0_flat = gamma0.divide(scf)\n",
        "        gamma0_flatDB = (ee.Image.constant(10)\n",
        "                         .multiply(gamma0_flat.log10())\n",
        "                         .select(['VV', 'VH'],['VV_gamma0flat', 'VH_gamma0flat'])\n",
        "                        )\n",
        "\n",
        "        masks = _masking(alpha_rRad, theta_iRad, buffer)\n",
        "\n",
        "        # calculate the ratio for RGB vis\n",
        "        ratio_flat = (gamma0_flatDB.select('VV_gamma0flat')\n",
        "                        .subtract(gamma0_flatDB.select('VH_gamma0flat'))\n",
        "                        .rename('ratio_gamma0flat')\n",
        "                     )\n",
        "\n",
        "        return (image.rename(['VV_sigma0', 'VH_sigma0', 'incAngle'])\n",
        "                      .addBands(gamma0dB)\n",
        "                      .addBands(ratio_gamma)\n",
        "                      .addBands(gamma0_flatDB)\n",
        "                      .addBands(ratio_flat)\n",
        "                      .addBands(alpha_rRad.rename('alpha_rRad'))\n",
        "                      .addBands(alpha_azRad.rename('alpha_azRad'))\n",
        "                      .addBands(phi_sRad.rename('aspect'))\n",
        "                      .addBands(alpha_sRad.rename('slope'))\n",
        "                      .addBands(theta_iRad.rename('theta_iRad'))\n",
        "                      .addBands(theta_liaRad.rename('theta_liaRad'))\n",
        "                      .addBands(masks)\n",
        "                      .addBands(height.rename('elevation'))\n",
        "                 )    \n",
        "    \n",
        "    # run and return correction\n",
        "    return collection.map(_correct)"
      ],
      "metadata": {
        "id": "KcJqUPrGfisY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "2. \n",
        "*   Defining the AOI\n",
        "*   Loading the UKECH 10m land cover file from Google Earth Engine (GEE) assets\n",
        "*   Filtering the Sentinel-1 imagery\n",
        "*   Loading the BlueSky 5m DEM from Google Earth Engine (GEE) assets\n",
        "*   Looping (Mapping) through the slope correction module calculations\n",
        "\n",
        "NOTE: The computations are done on the Google Server and take several hours to complete. At the end the computations are stored within respective folders created in the user's Google Drive \n",
        "\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "mEmMp3y6quxB"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# geometry for AOI in Hollin Hill North Yorkshire\n",
        "geometry = ee.Geometry.Polygon([[[-1.231403738843886, 54.291069000328001 ], [ -0.687709680649565, 54.289328250942106 ], [ -0.679296058617742, 53.914486883179535 ], [ -1.228212364969746, 53.913906633384236 ], [ -1.231403738843886, 54.291069000328001]]])\n",
        "\n",
        "# Corine land cover (2018) of Hollin HIll\n",
        "lc = ee.Image('users/kennethcassar/LC/LCM_10m').rename('landcover')\n",
        "\n",
        "# filter Sentinel-1 collection for AOI and selected dates\n",
        "s1Collection = ee.ImageCollection('COPERNICUS/S1_GRD') \\\n",
        "        .filterBounds(geometry) \\\n",
        "        .filterDate('2019-01-01', '2021-12-31')\n",
        "\n",
        "\n",
        "# path to dem \n",
        "dem = ee.Image('users/kennethcassar/DEM/Bluesky_DTM_Clip')\n",
        "\n",
        "# list of models\n",
        "models = ['volume', 'surface']\n",
        "\n",
        "# this is the scale we want the data to be sampled\n",
        "scale = 10 \n",
        "\n",
        "# loop through all combinations and export to drive\n",
        "for model in models:\n",
        "    \n",
        "    # get the respective collection and bands and mosaic to a single image\n",
        "    corrected_image = slope_correction(\n",
        "        s1Collection, \n",
        "        ee.Image(dem), \n",
        "        model\n",
        "    ).mosaic()\n",
        "\n",
        "    # we get geometry and projection of the image\n",
        "    proj = corrected_image.select(1).projection()\n",
        "    geom = corrected_image.clip(geometry).select(1).geometry()\n",
        "\n",
        "    # add lc and bring everything to same projection/geometry\n",
        "    added_LC = corrected_image.addBands(lc)\n",
        "    image_reprojected = added_LC.reproject(proj, scale=scale).clip(geom)\n",
        "      \n",
        "    # get the bandlist \n",
        "    bandlist = image_reprojected.bandNames().getInfo()\n",
        "    \n",
        "    # create an export job for each band\n",
        "    for band in bandlist:\n",
        "        task = ee.batch.Export.image.toDrive(\n",
        "            image=image_reprojected.select(band).clip(geom),\n",
        "            description=band,\n",
        "            folder='slope_correction_LCM/{}_{}_buf_0'.format(dem, model),\n",
        "            fileNamePrefix=band,\n",
        "            region=geom.coordinates().getInfo(),\n",
        "            scale=scale,\n",
        "            maxPixels=1e12\n",
        "            )\n",
        "        task.start()\n",
        "        print(task.status())"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nif85tcnu3Pt",
        "outputId": "2d3ae87e-e3f2-4f67-d660-3c5a47ce6aec"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "{'state': 'READY', 'description': 'VV_sigma0', 'creation_timestamp_ms': 1657312678121, 'update_timestamp_ms': 1657312678121, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'EUJE4VUIJIKU524LZBLDE4I4', 'name': 'projects/earthengine-legacy/operations/EUJE4VUIJIKU524LZBLDE4I4'}\n",
            "{'state': 'READY', 'description': 'VH_sigma0', 'creation_timestamp_ms': 1657312679164, 'update_timestamp_ms': 1657312679164, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'G4ZLBVABWZUBUSF4426YA45V', 'name': 'projects/earthengine-legacy/operations/G4ZLBVABWZUBUSF4426YA45V'}\n",
            "{'state': 'READY', 'description': 'incAngle', 'creation_timestamp_ms': 1657312680230, 'update_timestamp_ms': 1657312680230, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'VZBRHQ22PUDBSNXXMB3RQ2QC', 'name': 'projects/earthengine-legacy/operations/VZBRHQ22PUDBSNXXMB3RQ2QC'}\n",
            "{'state': 'READY', 'description': 'VV_gamma0', 'creation_timestamp_ms': 1657312681020, 'update_timestamp_ms': 1657312681020, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'PTTXDU3AXBSSR4IIAAXNJX3K', 'name': 'projects/earthengine-legacy/operations/PTTXDU3AXBSSR4IIAAXNJX3K'}\n",
            "{'state': 'READY', 'description': 'VH_gamma0', 'creation_timestamp_ms': 1657312681806, 'update_timestamp_ms': 1657312681806, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'DTUJAFWNJCXV3JMXBELMS2O6', 'name': 'projects/earthengine-legacy/operations/DTUJAFWNJCXV3JMXBELMS2O6'}\n",
            "{'state': 'READY', 'description': 'ratio_gamma0', 'creation_timestamp_ms': 1657312682630, 'update_timestamp_ms': 1657312682630, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'ODT5PUU62VUVX7JQRTG3AHU6', 'name': 'projects/earthengine-legacy/operations/ODT5PUU62VUVX7JQRTG3AHU6'}\n",
            "{'state': 'READY', 'description': 'VV_gamma0flat', 'creation_timestamp_ms': 1657312683456, 'update_timestamp_ms': 1657312683456, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'TBXDY34GCXCMBU4OJBZ55S2L', 'name': 'projects/earthengine-legacy/operations/TBXDY34GCXCMBU4OJBZ55S2L'}\n",
            "{'state': 'READY', 'description': 'VH_gamma0flat', 'creation_timestamp_ms': 1657312684260, 'update_timestamp_ms': 1657312684260, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'VW6IQWVBRTJZWXHF2CU7GUI4', 'name': 'projects/earthengine-legacy/operations/VW6IQWVBRTJZWXHF2CU7GUI4'}\n",
            "{'state': 'READY', 'description': 'ratio_gamma0flat', 'creation_timestamp_ms': 1657312685119, 'update_timestamp_ms': 1657312685119, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'SYIBKL74SF4OKIMNW4NJKTTJ', 'name': 'projects/earthengine-legacy/operations/SYIBKL74SF4OKIMNW4NJKTTJ'}\n",
            "{'state': 'READY', 'description': 'alpha_rRad', 'creation_timestamp_ms': 1657312686015, 'update_timestamp_ms': 1657312686015, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'M77IQCD6MQG7MOPI5YEMXIPA', 'name': 'projects/earthengine-legacy/operations/M77IQCD6MQG7MOPI5YEMXIPA'}\n",
            "{'state': 'READY', 'description': 'alpha_azRad', 'creation_timestamp_ms': 1657312686853, 'update_timestamp_ms': 1657312686853, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'TTGFECXBMB4BE7HHBEKGM6XY', 'name': 'projects/earthengine-legacy/operations/TTGFECXBMB4BE7HHBEKGM6XY'}\n",
            "{'state': 'READY', 'description': 'aspect', 'creation_timestamp_ms': 1657312687705, 'update_timestamp_ms': 1657312687705, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'HOO2N6LBLOJAOGMDL3A7ACOF', 'name': 'projects/earthengine-legacy/operations/HOO2N6LBLOJAOGMDL3A7ACOF'}\n",
            "{'state': 'READY', 'description': 'slope', 'creation_timestamp_ms': 1657312688552, 'update_timestamp_ms': 1657312688552, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': '42YLMMP2EMSAKHIBDKT6SMMD', 'name': 'projects/earthengine-legacy/operations/42YLMMP2EMSAKHIBDKT6SMMD'}\n",
            "{'state': 'READY', 'description': 'theta_iRad', 'creation_timestamp_ms': 1657312689397, 'update_timestamp_ms': 1657312689397, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'OXAGNAV5GW4E5GTBWDQTVQDK', 'name': 'projects/earthengine-legacy/operations/OXAGNAV5GW4E5GTBWDQTVQDK'}\n",
            "{'state': 'READY', 'description': 'theta_liaRad', 'creation_timestamp_ms': 1657312690226, 'update_timestamp_ms': 1657312690226, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'ADWRSB5ZR76AOSAJFGRF6PVA', 'name': 'projects/earthengine-legacy/operations/ADWRSB5ZR76AOSAJFGRF6PVA'}\n",
            "{'state': 'READY', 'description': 'layover', 'creation_timestamp_ms': 1657312690973, 'update_timestamp_ms': 1657312690973, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': '4WURCHUA26YLAYTUXGA7353S', 'name': 'projects/earthengine-legacy/operations/4WURCHUA26YLAYTUXGA7353S'}\n",
            "{'state': 'READY', 'description': 'shadow', 'creation_timestamp_ms': 1657312691756, 'update_timestamp_ms': 1657312691756, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'AVMJDIOPTZWKGHJYQSUMH7LK', 'name': 'projects/earthengine-legacy/operations/AVMJDIOPTZWKGHJYQSUMH7LK'}\n",
            "{'state': 'READY', 'description': 'no_data_mask', 'creation_timestamp_ms': 1657312692707, 'update_timestamp_ms': 1657312692707, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'XQI6NMMI7XODASJ55QSCGHK7', 'name': 'projects/earthengine-legacy/operations/XQI6NMMI7XODASJ55QSCGHK7'}\n",
            "{'state': 'READY', 'description': 'elevation', 'creation_timestamp_ms': 1657312693553, 'update_timestamp_ms': 1657312693553, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': '2EQ5LRAIKIPOWQAPBS4CWEJU', 'name': 'projects/earthengine-legacy/operations/2EQ5LRAIKIPOWQAPBS4CWEJU'}\n",
            "{'state': 'READY', 'description': 'landcover', 'creation_timestamp_ms': 1657312694420, 'update_timestamp_ms': 1657312694420, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'BYTMXYRUM3N4MX7C3S3IZSAG', 'name': 'projects/earthengine-legacy/operations/BYTMXYRUM3N4MX7C3S3IZSAG'}\n",
            "{'state': 'READY', 'description': 'VV_sigma0', 'creation_timestamp_ms': 1657312695614, 'update_timestamp_ms': 1657312695614, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'TFZ6VCPAIH3OQTQFFINW5NBC', 'name': 'projects/earthengine-legacy/operations/TFZ6VCPAIH3OQTQFFINW5NBC'}\n",
            "{'state': 'READY', 'description': 'VH_sigma0', 'creation_timestamp_ms': 1657312696407, 'update_timestamp_ms': 1657312696407, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'XEHR4DTCX3QO45WXQWLX7MDJ', 'name': 'projects/earthengine-legacy/operations/XEHR4DTCX3QO45WXQWLX7MDJ'}\n",
            "{'state': 'READY', 'description': 'incAngle', 'creation_timestamp_ms': 1657312697183, 'update_timestamp_ms': 1657312697183, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'GONRCYQBENX5KQBDAXW5IJCD', 'name': 'projects/earthengine-legacy/operations/GONRCYQBENX5KQBDAXW5IJCD'}\n",
            "{'state': 'READY', 'description': 'VV_gamma0', 'creation_timestamp_ms': 1657312698616, 'update_timestamp_ms': 1657312698616, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'WQJNVJKMWMHQIT3EQIQA3RF3', 'name': 'projects/earthengine-legacy/operations/WQJNVJKMWMHQIT3EQIQA3RF3'}\n",
            "{'state': 'READY', 'description': 'VH_gamma0', 'creation_timestamp_ms': 1657312699450, 'update_timestamp_ms': 1657312699450, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'UNS52ONG7VU3KBJYFLM3Y257', 'name': 'projects/earthengine-legacy/operations/UNS52ONG7VU3KBJYFLM3Y257'}\n",
            "{'state': 'READY', 'description': 'ratio_gamma0', 'creation_timestamp_ms': 1657312700286, 'update_timestamp_ms': 1657312700286, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'PLNNEV5BIPQO5W2PHZFYEG3T', 'name': 'projects/earthengine-legacy/operations/PLNNEV5BIPQO5W2PHZFYEG3T'}\n",
            "{'state': 'READY', 'description': 'VV_gamma0flat', 'creation_timestamp_ms': 1657312701051, 'update_timestamp_ms': 1657312701051, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'DQMHFKVHZ6X4SBY5NZSOG6UH', 'name': 'projects/earthengine-legacy/operations/DQMHFKVHZ6X4SBY5NZSOG6UH'}\n",
            "{'state': 'READY', 'description': 'VH_gamma0flat', 'creation_timestamp_ms': 1657312701857, 'update_timestamp_ms': 1657312701857, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': '54MBLFXVXOJRKDTGAKBJHLZW', 'name': 'projects/earthengine-legacy/operations/54MBLFXVXOJRKDTGAKBJHLZW'}\n",
            "{'state': 'READY', 'description': 'ratio_gamma0flat', 'creation_timestamp_ms': 1657312702723, 'update_timestamp_ms': 1657312702723, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'AW2LXXBZQJCHMHI5JP6UBFOM', 'name': 'projects/earthengine-legacy/operations/AW2LXXBZQJCHMHI5JP6UBFOM'}\n",
            "{'state': 'READY', 'description': 'alpha_rRad', 'creation_timestamp_ms': 1657312703625, 'update_timestamp_ms': 1657312703625, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'CVFKTHPWW6Q7WXXFSXJXZ6D6', 'name': 'projects/earthengine-legacy/operations/CVFKTHPWW6Q7WXXFSXJXZ6D6'}\n",
            "{'state': 'READY', 'description': 'alpha_azRad', 'creation_timestamp_ms': 1657312704453, 'update_timestamp_ms': 1657312704453, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'CCDLXW4E74ICWPR7X7SJDF3E', 'name': 'projects/earthengine-legacy/operations/CCDLXW4E74ICWPR7X7SJDF3E'}\n",
            "{'state': 'READY', 'description': 'aspect', 'creation_timestamp_ms': 1657312705212, 'update_timestamp_ms': 1657312705212, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'GWOJLCMQKR6SRZZ3PE36AGZ6', 'name': 'projects/earthengine-legacy/operations/GWOJLCMQKR6SRZZ3PE36AGZ6'}\n",
            "{'state': 'READY', 'description': 'slope', 'creation_timestamp_ms': 1657312706064, 'update_timestamp_ms': 1657312706064, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'ARJEUYUWHG2C3CYL4F3PYKGP', 'name': 'projects/earthengine-legacy/operations/ARJEUYUWHG2C3CYL4F3PYKGP'}\n",
            "{'state': 'READY', 'description': 'theta_iRad', 'creation_timestamp_ms': 1657312706865, 'update_timestamp_ms': 1657312706865, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'OYKBXQWGT2HTPMUC7ETSCOJT', 'name': 'projects/earthengine-legacy/operations/OYKBXQWGT2HTPMUC7ETSCOJT'}\n",
            "{'state': 'READY', 'description': 'theta_liaRad', 'creation_timestamp_ms': 1657312707753, 'update_timestamp_ms': 1657312707753, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'ZUHGLMF53ZCZZ2LGB3GPP4MX', 'name': 'projects/earthengine-legacy/operations/ZUHGLMF53ZCZZ2LGB3GPP4MX'}\n",
            "{'state': 'READY', 'description': 'layover', 'creation_timestamp_ms': 1657312708542, 'update_timestamp_ms': 1657312708542, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'F4AWGREADHAKXBAKZS7JH3DL', 'name': 'projects/earthengine-legacy/operations/F4AWGREADHAKXBAKZS7JH3DL'}\n",
            "{'state': 'READY', 'description': 'shadow', 'creation_timestamp_ms': 1657312709290, 'update_timestamp_ms': 1657312709290, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'M3TC46UQQVZCLCKZLO5W5GYG', 'name': 'projects/earthengine-legacy/operations/M3TC46UQQVZCLCKZLO5W5GYG'}\n",
            "{'state': 'READY', 'description': 'no_data_mask', 'creation_timestamp_ms': 1657312710087, 'update_timestamp_ms': 1657312710087, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'SQCDD4SWDUQK347FCW5X5FDT', 'name': 'projects/earthengine-legacy/operations/SQCDD4SWDUQK347FCW5X5FDT'}\n",
            "{'state': 'READY', 'description': 'elevation', 'creation_timestamp_ms': 1657312710954, 'update_timestamp_ms': 1657312710954, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'FZYKGNW6SMKEFOGAFQ44PXR7', 'name': 'projects/earthengine-legacy/operations/FZYKGNW6SMKEFOGAFQ44PXR7'}\n",
            "{'state': 'READY', 'description': 'landcover', 'creation_timestamp_ms': 1657312711833, 'update_timestamp_ms': 1657312711833, 'start_timestamp_ms': 0, 'task_type': 'EXPORT_IMAGE', 'id': 'ZPH4IBCMVHNE7XIGCSQWKMVM', 'name': 'projects/earthengine-legacy/operations/ZPH4IBCMVHNE7XIGCSQWKMVM'}\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "3. Downloading of necessary Python libraries"
      ],
      "metadata": {
        "id": "5ajZXciMrXYA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!apt-get -qq install -y libspatialindex-dev\n",
        "!pip install -q --upgrade earthpy"
      ],
      "metadata": {
        "id": "Cmou-3-VsHK_",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "bc545086-c388-4bbb-8e22-c8b29abae09b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Selecting previously unselected package libspatialindex4v5:amd64.\n",
            "(Reading database ... 155639 files and directories currently installed.)\n",
            "Preparing to unpack .../libspatialindex4v5_1.8.5-5_amd64.deb ...\n",
            "Unpacking libspatialindex4v5:amd64 (1.8.5-5) ...\n",
            "Selecting previously unselected package libspatialindex-c4v5:amd64.\n",
            "Preparing to unpack .../libspatialindex-c4v5_1.8.5-5_amd64.deb ...\n",
            "Unpacking libspatialindex-c4v5:amd64 (1.8.5-5) ...\n",
            "Selecting previously unselected package libspatialindex-dev:amd64.\n",
            "Preparing to unpack .../libspatialindex-dev_1.8.5-5_amd64.deb ...\n",
            "Unpacking libspatialindex-dev:amd64 (1.8.5-5) ...\n",
            "Setting up libspatialindex4v5:amd64 (1.8.5-5) ...\n",
            "Setting up libspatialindex-c4v5:amd64 (1.8.5-5) ...\n",
            "Setting up libspatialindex-dev:amd64 (1.8.5-5) ...\n",
            "Processing triggers for libc-bin (2.27-3ubuntu1.3) ...\n",
            "/sbin/ldconfig.real: /usr/local/lib/python3.7/dist-packages/ideep4py/lib/libmkldnn.so.0 is not a symbolic link\n",
            "\n",
            "\u001b[K     |████████████████████████████████| 1.4 MB 4.3 MB/s \n",
            "\u001b[K     |████████████████████████████████| 1.0 MB 36.9 MB/s \n",
            "\u001b[K     |████████████████████████████████| 19.3 MB 40.1 MB/s \n",
            "\u001b[K     |████████████████████████████████| 6.3 MB 16.0 MB/s \n",
            "\u001b[K     |████████████████████████████████| 16.7 MB 15.9 MB/s \n",
            "\u001b[?25h"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "4. Access the Google Drive path to eventually save the results in "
      ],
      "metadata": {
        "id": "i1UFDTmYrfhJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "metadata": {
        "id": "rFsVg1phsSHo",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "a85a4166-f74a-49d5-8cf1-92c58cdbc579"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "import itertools\n",
        "\n",
        "import numpy as np\n",
        "import pandas as pd \n",
        "from scipy import stats\n",
        "from scipy import optimize\n",
        "\n",
        "import rasterio as rio\n",
        "from rasterio.windows import Window\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import earthpy.plot as ep\n",
        "%matplotlib inline"
      ],
      "metadata": {
        "id": "XgiEpuO4sXKG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "5. Read the necessary bands calculated under step 2"
      ],
      "metadata": {
        "id": "ZZM8-6s-ryMh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "list_of_layers = [\n",
        "    'VV_gamma0', 'VH_gamma0',\n",
        "    'VV_gamma0flat', 'VH_gamma0flat',\n",
        "    'alpha_rRad', 'theta_liaRad', 'aspect', \n",
        "    'layover', 'shadow',     \n",
        "    'landcover']\n",
        "\n",
        "# paths to dem and model types\n",
        "dem = 'Bluesky_DTM_Clip'\n",
        "models = ['volume', 'surface']\n",
        "buffers = [0]\n",
        "\n",
        "modelDict = {}\n",
        "\n",
        "# loop thorugh all combinations and put into dictionary\n",
        "for model, buffer in itertools.product(models, buffers):\n",
        "\n",
        "    key = '{}_{}_buf_{}'.format(dem, model, buffer)\n",
        "    # here we read all layers into a dictionary\n",
        "    dataDict = {}\n",
        "    for layer in list_of_layers:\n",
        "        with rio.open('/content/drive/My Drive/slope_correction_LCM {}/{}.tif'.format(key, layer)) as src:\n",
        "            print('Reading ' + layer)\n",
        "            dataDict[layer] = np.nan_to_num(src.read(window=Window(0, 340, 3000, 3000)))[0]\n",
        "            print(dataDict[layer].shape)\n",
        "    \n",
        "    # write respective dataDict to our model dict, where different models are stored\n",
        "    modelDict[key] = dataDict"
      ],
      "metadata": {
        "id": "62HXEE9isZ4e",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c5213e9e-1be0-457b-86dc-6c6ff2e91e2d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Reading VV_gamma0\n",
            "(3000, 3000)\n",
            "Reading VH_gamma0\n",
            "(3000, 3000)\n",
            "Reading VV_gamma0flat\n",
            "(3000, 3000)\n",
            "Reading VH_gamma0flat\n",
            "(3000, 3000)\n",
            "Reading alpha_rRad\n",
            "(3000, 3000)\n",
            "Reading theta_liaRad\n",
            "(3000, 3000)\n",
            "Reading aspect\n",
            "(3000, 3000)\n",
            "Reading layover\n",
            "(3000, 3000)\n",
            "Reading shadow\n",
            "(3000, 3000)\n",
            "Reading landcover\n",
            "(3000, 3000)\n",
            "Reading VV_gamma0\n",
            "(3000, 3000)\n",
            "Reading VH_gamma0\n",
            "(3000, 3000)\n",
            "Reading VV_gamma0flat\n",
            "(3000, 3000)\n",
            "Reading VH_gamma0flat\n",
            "(3000, 3000)\n",
            "Reading alpha_rRad\n",
            "(3000, 3000)\n",
            "Reading theta_liaRad\n",
            "(3000, 3000)\n",
            "Reading aspect\n",
            "(3000, 3000)\n",
            "Reading layover\n",
            "(3000, 3000)\n",
            "Reading shadow\n",
            "(3000, 3000)\n",
            "Reading landcover\n",
            "(3000, 3000)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "6. Calculate the statistics and plot the results as charts. At this stage the charts are just stored as temporary files and not visualised."
      ],
      "metadata": {
        "id": "Sg2ToIEEsKlU"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def create_plot_aspect_against_backscatter(model, data, array, mask, stats_dict, outfile, gridsize):\n",
        "    \n",
        "    # getlayer info #\n",
        "    polarisation = layer.split('_')[0]\n",
        "    flat = layer.split('_')[1][-4:]\n",
        "    if flat != 'flat':\n",
        "        flat = False\n",
        "    \n",
        "    # fore and backslope lines\n",
        "    look_angle = 33.089022518635176 # hardcoded \n",
        "    backslope = look_angle\n",
        "    foreslope = backslope + 180\n",
        "    vertical_y = np.linspace(-28, 8, 5)\n",
        "    fs_x = 0 * vertical_y + foreslope\n",
        "    bs_x = 0 * vertical_y + backslope\n",
        "    \n",
        "    # calculate mean line\n",
        "    horizontal_x = np.linspace(0, 360, 10) \n",
        "    mean_y = 0 * horizontal_x + stats_dict['mean']\n",
        "    \n",
        "    if not flat:\n",
        "        y_label = r'$\\gamma^0$ [dB]'\n",
        "    else:\n",
        "        y_label = r'$\\gamma^0_f$ [dB]'\n",
        "    \n",
        "    # check for 0s in aspect and update mask\n",
        "    data['aspect'][data['aspect'] == 0] = np.nan\n",
        "    mask = mask & np.isfinite(data['aspect'])\n",
        "    aspect_deg_masked = np.subtract(to_deg(data['aspect'][mask]), 180)\n",
        "    aspect_deg_masked = to_deg(data['aspect'][mask])\n",
        "\n",
        "    # plot\n",
        "    # surpress plotting, since we only want to save the files\n",
        "    plt.ioff()\n",
        "    X = sns.jointplot(aspect_deg_masked, array[mask], kind='hex', gridsize=(gridsize, gridsize))\n",
        "    X.ax_joint.plot(horizontal_x, mean_y, 'k--', linewidth=.75)\n",
        "    X.ax_joint.plot(fs_x, vertical_y, 'k--', linewidth=.75)\n",
        "    X.ax_joint.plot(bs_x, vertical_y, 'k--', linewidth=.75)\n",
        "    X.ax_joint.set_xlabel(r'Aspect angle $\\phi_s$ [deg]', fontsize=14)\n",
        "    X.ax_joint.set_ylabel(y_label,  fontsize=14)\n",
        "    X.ax_joint.set_ylim(-30, 10)\n",
        "    #X.ax_joint.set_xlim(-190, 185)\n",
        "    X.ax_joint.set_xlim(-10, 365)\n",
        "    \n",
        "    # add textbox with ampl, mean and sd\n",
        "    props = dict(boxstyle='round', facecolor='lightgrey', alpha=0.5)\n",
        "    textstr = '\\n'.join((\n",
        "        r'$\\mathrm{A}=%.2f$' % (stats_dict['amplitude'], ),\n",
        "        r'$\\mu=%.2f$' % (stats_dict['mean'], ),\n",
        "        r'$\\sigma=%.2f$' % (stats_dict['sd'], )))\n",
        "    X.ax_joint.text(290, -28,textstr, fontsize=10, bbox=props)\n",
        "    \n",
        "    # add title\n",
        "    if not flat:\n",
        "        title = '{} Backscatter modulation by slopes'.format(polarisation)\n",
        "    else:\n",
        "        title = '{} Backscatter after Model {}'.format(polarisation, model)\n",
        "        \n",
        "    plt.suptitle(title, x=0.45, y=1.02, fontweight='bold', fontsize=14)\n",
        "    \n",
        "    # save\n",
        "    plt.savefig(outfile, bbox_inches='tight', pad_inches=0.5)\n",
        "    plt.close()"
      ],
      "metadata": {
        "id": "9VklpDg4siTP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def create_plot_alpha_against_backscatter(model, data, array, mask, stats_dict, outfile, gridsize):\n",
        "    \n",
        "    # getlayer info #\n",
        "    polarisation = layer.split('_')[0]\n",
        "    flat = layer.split('_')[1][-4:]\n",
        "    if flat != 'flat':\n",
        "        flat = False\n",
        "    \n",
        "    # calculate mean line\n",
        "    horizontal_x = np.linspace(-40, 40, 10) \n",
        "    mean_y = 0 * horizontal_x + stats_dict['mean']\n",
        "    \n",
        "    if not flat:\n",
        "        y_label = r'$\\gamma^0$ [dB]'\n",
        "    else:\n",
        "        y_label = r'$\\gamma^0_f$ [dB]'\n",
        "    \n",
        "    # check for 0s in aspect and update mask\n",
        "    data['alpha_rRad'][data['alpha_rRad'] == 0] = np.nan\n",
        "    mask = mask & np.isfinite(data['alpha_rRad'])\n",
        "    alpha_deg_masked = to_deg(data['alpha_rRad'][mask])\n",
        "\n",
        "    # plot\n",
        "    # surpress plotting, since we only want to save the files\n",
        "    plt.ioff()\n",
        "    X = sns.jointplot(alpha_deg_masked, array[mask], kind='hex', gridsize=(gridsize, gridsize))\n",
        "    X.ax_joint.plot(horizontal_x, mean_y, 'k--', linewidth=.75)\n",
        "    X.ax_joint.set_xlabel(r'Slope Steepness in range $\\alpha$ [deg]', fontsize=14)\n",
        "    X.ax_joint.set_ylabel(y_label, fontsize=14)\n",
        "    X.ax_joint.set_ylim(-30, 10)\n",
        "    X.ax_joint.set_xlim(-45, 45)\n",
        "    \n",
        "    # add textbox with ampl, mean and sd\n",
        "    props = dict(boxstyle='round', facecolor='lightgrey', alpha=0.5)\n",
        "    textstr = '\\n'.join((\n",
        "        r'$\\mathrm{s}=%.2f$' % (stats_dict['slope'], ),\n",
        "        r'$\\mu=%.2f$' % (stats_dict['mean'], ),\n",
        "        r'$\\sigma=%.2f$' % (stats_dict['sd'], )))\n",
        "    X.ax_joint.text(20, -28,textstr, fontsize=10, bbox=props)\n",
        "    \n",
        "    # add title\n",
        "    if not flat:\n",
        "        title = '{} Backscatter modulation by slopes'.format(polarisation)\n",
        "    else:\n",
        "        title = '{} Backscatter after Model {}'.format(polarisation, model)\n",
        "        \n",
        "    plt.suptitle(title, x=0.45, y=1.02, fontweight='bold', fontsize=14)\n",
        "    \n",
        "    # save\n",
        "    plt.savefig(outfile, bbox_inches='tight', pad_inches=0.5)\n",
        "    plt.close()"
      ],
      "metadata": {
        "id": "rtm_OWjqsmpb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def create_plot_lia_against_backscatter(model, data, array, mask, stats_dict, outfile, gridsize):\n",
        "    \n",
        "    \n",
        "    # getlayer info #\n",
        "    polarisation = layer.split('_')[0]\n",
        "    flat = layer.split('_')[1][-4:]\n",
        "    if flat != 'flat':\n",
        "        flat = False\n",
        "    \n",
        "    # calculate mean line\n",
        "    horizontal_x = np.linspace(0, 90, 10) \n",
        "    mean_y = 0 * horizontal_x + stats_dict['mean']\n",
        "    \n",
        "    if not flat:\n",
        "        y_label = r'$\\gamma^0$ [dB]'\n",
        "    else:\n",
        "        y_label = r'$\\gamma^0_f$ [dB]'\n",
        "    \n",
        "    # check for 0s in aspect and update mask\n",
        "    data['theta_liaRad'][data['theta_liaRad'] == 0] = np.nan\n",
        "    mask = mask & np.isfinite(data['theta_liaRad'])\n",
        "    theta_deg_masked = to_deg(data['theta_liaRad'][mask])\n",
        "\n",
        "    # plot\n",
        "    # surpress plotting, since we only want to save the files\n",
        "    plt.ioff()\n",
        "    X = sns.jointplot(theta_deg_masked, array[mask], kind='hex', gridsize=(gridsize, gridsize))\n",
        "    X.ax_joint.plot(horizontal_x, mean_y, 'k--', linewidth=.75)\n",
        "    X.ax_joint.set_xlabel(r'Local Incidence Angle $\\theta$ [deg]', fontsize=14)\n",
        "    X.ax_joint.set_ylabel(y_label, fontsize=14)\n",
        "    X.ax_joint.set_ylim(-30, 10)\n",
        "    X.ax_joint.set_xlim(-10, 100)\n",
        "    \n",
        "    # add textbox with ampl, mean and sd\n",
        "    props = dict(boxstyle='round', facecolor='lightgrey', alpha=0.5)\n",
        "    textstr = '\\n'.join((\n",
        "        #r'$\\mathrm{s}=%.2f$' % (stats_dict['slope'], ),\n",
        "        r'$\\mu=%.2f$' % (stats_dict['mean'], ),\n",
        "        r'$\\sigma=%.2f$' % (stats_dict['sd'], )))\n",
        "    X.ax_joint.text(-5, -28,textstr, fontsize=10, bbox=props)\n",
        "    \n",
        "    # add title\n",
        "    if not flat:\n",
        "        title = '{} Backscatter modulation by slopes'.format(polarisation)\n",
        "    else:\n",
        "        title = '{} Backscatter after Model {}'.format(polarisation, model)\n",
        "        \n",
        "    plt.suptitle(title, x=0.45, y=1.02, fontweight='bold', fontsize=14)\n",
        "    \n",
        "    # save\n",
        "    plt.savefig(outfile, bbox_inches='tight', pad_inches=0.5)\n",
        "    plt.close()"
      ],
      "metadata": {
        "id": "NeIX4ha9sqpU"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# from rad to degree\n",
        "def to_deg(rad):\n",
        "    return np.multiply(rad,  np.divide(180,np.pi))\n",
        "\n",
        "# sin function or curve fitting\n",
        "def sin_func(x, a, b, c):\n",
        "    return a * np.sin(x + b) - c\n",
        "\n",
        "def tf_stats(data, array, layer, lc_class, mask):\n",
        "    \n",
        "    #------------------------------------------------\n",
        "    # slope calculation\n",
        "    \n",
        "    # mask alpha angle\n",
        "    data['alpha_rRad'][data['alpha_rRad'] == 0] = np.nan\n",
        "    mask_alpha = mask & np.isfinite(data['alpha_rRad'])\n",
        "    alpha_deg_masked = np.subtract(to_deg(data['alpha_rRad'][mask_alpha]), 180)\n",
        "    \n",
        "    # mask out nans\n",
        "    x = array.copy()\n",
        "    x[~mask_alpha] = np.nan\n",
        "    x = x[np.logical_not(np.isnan(x))] \n",
        "    y = alpha_deg_masked\n",
        "    y = y[np.logical_not(np.isnan(y))]\n",
        "    \n",
        "    # lin-regression\n",
        "    slope, intercept, r_value, p_value, std_err = stats.linregress(y, x)        \n",
        "    #------------------------------------------------\n",
        "    \n",
        "    #------------------------------------------------\n",
        "    # amplitude calculation\n",
        "    \n",
        "    # mask aspect 0s and nans\n",
        "    data['aspect'][data['aspect'] == 0] = np.nan\n",
        "    mask_aspect = mask & np.isfinite(data['aspect'])\n",
        "    aspect_deg_masked = np.subtract(data['aspect'][mask_aspect], np.pi)\n",
        "    \n",
        "    # mask out nans\n",
        "    x = array.copy()\n",
        "    x[~mask_aspect] = np.nan\n",
        "    x = x[np.logical_not(np.isnan(x))]\n",
        "    y = aspect_deg_masked\n",
        "    y = y[np.logical_not(np.isnan(y))]\n",
        "   \n",
        "    # curve fitting\n",
        "    params, params_covariance = optimize.curve_fit(sin_func, y, x)\n",
        "    amp = params[0]\n",
        "    #------------------------------------------------\n",
        "\n",
        "    # mean, sd\n",
        "    mean = np.nanmean(array[mask_alpha])\n",
        "    std = np.nanstd(array[mask_alpha])\n",
        "    \n",
        "    # create final dictionary\n",
        "    stat_dict = {'lc_class': lc_class, \n",
        "                 'layer': layer, \n",
        "                 'count': np.sum(mask), \n",
        "                 'mean': mean, \n",
        "                 'sd': std, \n",
        "                 'slope': slope, \n",
        "                 'amplitude': np.abs(amp)\n",
        "                }\n",
        "    \n",
        "    return stat_dict"
      ],
      "metadata": {
        "id": "JdHDkZg2sv3G"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "7. Define the land cover classes according to the loaded land cover file under step 2. The calculated statistics are then mapped against the land cover classes."
      ],
      "metadata": {
        "id": "UEd5QHvDssFw"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# list of class names\n",
        "legend_entries = [\"Broad-leaved forest\", \"Arable & Horticulture\", \"Improved Grassland\", \"Calcareous Grassland\", \"Heather Grassland\"]\n",
        "# list of class values\n",
        "ras_values = [1, 3, 4, 6, 10]\n",
        "\n",
        "# prepare columns for new dataframe\n",
        "df_cols = ['lc_class', 'layer', 'count', 'mean', 'sd', 'slope', 'amplitude']\n",
        "\n",
        "# paths to dem and model types and buffer\n",
        "dem = 'Bluesky_DTM_Clip'\n",
        "models = ['surface', 'volume']\n",
        "buffers = [0] \n",
        "\n",
        "# loop through different models and buffers\n",
        "for model, buffer in itertools.product(models, buffers):\n",
        "    \n",
        "    # get respective arrays within the model/datadict\n",
        "    key = '{}_{}_buf_{}'.format(dem, model, buffer)\n",
        "    dataDict = modelDict[key]\n",
        "    \n",
        "    print(' INFO: Creating figures and stats for {} {} with buffer {}.'\n",
        "      .format(dem, model, buffer)\n",
        "    )\n",
        "    # create empty dataframe for statistics\n",
        "    df_stats = pd.DataFrame(columns=df_cols)\n",
        "    \n",
        "    # crate outdirectory where plots and stats will be saved\n",
        "    outdir = '/content/drive/My Drive/slope_correction_LCM/pictures/{}/'.format(key)\n",
        "    os.makedirs(outdir, exist_ok=True)\n",
        "    \n",
        "    # loop through classes and respective raster values file\n",
        "    for legend_entry, ras_value in zip(legend_entries, ras_values):\n",
        "        \n",
        "        # set raster value respective to class\n",
        "        print(' INFO: Analysing {} with raster value {}.'\n",
        "          .format(legend_entry, ras_value)\n",
        "        )\n",
        "\n",
        "        # loop thorugh different corrected and uncorrected layers\n",
        "        for layer in ['VV_gamma0', 'VV_gamma0flat', 'VH_gamma0', 'VH_gamma0flat']:\n",
        "                        \n",
        "            # create combined Land Cover and Layover/Shadow mask\n",
        "            valid_data_mask = (\n",
        "                [dataDict['landcover'] == ras_value] & \n",
        "                (dataDict['layover'] > 0) & \n",
        "                (dataDict['shadow'] > 0)\n",
        "            )[0] \n",
        "\n",
        "            # apply this mask and add valid data mask of backscatter array \n",
        "            array = dataDict[layer].copy()\n",
        "            array[array == 0] = np.nan\n",
        "            mask = valid_data_mask & np.isfinite(array)\n",
        "\n",
        "            # set everything else to nan\n",
        "            array[~mask] = np.nan\n",
        "\n",
        "            # for some classes array might be empty, so we add an if\n",
        "            if True in np.unique(np.isfinite(array)):\n",
        "                \n",
        "                # calculate stats\n",
        "                stats_dict = tf_stats(\n",
        "                    dataDict.copy(), array, layer, legend_entry, mask\n",
        "                )\n",
        "                stats_dict['lc_class_code'] =  ras_value\n",
        "\n",
        "                # and put into pandas dataframe\n",
        "                df = pd.DataFrame([stats_dict], columns=stats_dict.keys())\n",
        "                df_stats = df_stats.append(stats_dict, ignore_index=True)\n",
        "\n",
        "                # plotting\n",
        "                gridsize=100\n",
        "                model_nr = '1' if model == 'volume' else '2'\n",
        "\n",
        "                create_plot_aspect_against_backscatter(\n",
        "                    model_nr, dataDict.copy(), array, mask, stats_dict, \n",
        "                    '{}/{}_{}_vs_aspect.png'.format(outdir, legend_entry, layer), \n",
        "                    gridsize\n",
        "                )\n",
        "\n",
        "                create_plot_alpha_against_backscatter(\n",
        "                    model_nr, dataDict.copy(), array, mask, stats_dict, \n",
        "                    '{}/{}_{}_vs_slope.png'.format(outdir, legend_entry, layer), \n",
        "                    gridsize\n",
        "                )\n",
        "\n",
        "                create_plot_lia_against_backscatter(\n",
        "                    model_nr, dataDict.copy(), array, mask, stats_dict, \n",
        "                    '{}/{}_{}_vs_LIA.png'.format(outdir, legend_entry, layer), \n",
        "                    gridsize\n",
        "                )\n",
        "\n",
        "        # save the complete stas dataframe to pickle\n",
        "        df_stats.reset_index()\n",
        "        df_stats.to_pickle('{}/stats.pickle'.format(outdir)) "
      ],
      "metadata": {
        "id": "b5i9Ud_fs0z9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "8. Save and visualise the statistical data computed for each land cover class within the AOI. "
      ],
      "metadata": {
        "id": "Lg5sG8J0tGI5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "df_stats_dict = {}\n",
        "\n",
        "# paths to dem and model types\n",
        "dem = 'Bluesky_DTM_Clip'\n",
        "models = ['surface', 'volume'] \n",
        "\n",
        "# this is for the concatenation of VV and VH stats\n",
        "concat_cols = ['lc_class', 'VV', 'Pixel count', 'VV mean', 'VV SD', 'VV slope', 'VV amp', \n",
        "                'VH', 'VH mean', 'VH SD', 'VH slope', 'VH amp']\n",
        "\n",
        "# create empty list for adding all stats within the for loop\n",
        "df_con_merged = []\n",
        "\n",
        "# loop through different dems and models\n",
        "for model in models:\n",
        "\n",
        "    # get respective arrays within the datadict\n",
        "    key = '{}_{}_buf_0'.format(dem, model)\n",
        "    \n",
        "    # read each df into the dictionary\n",
        "    outdir = '/content/drive/My Drive/slope_correction_LCM/pictures/{}/'.format(key)\n",
        "    df_stats = pd.read_pickle('{}/stats.pickle'.format(outdir))\n",
        "    \n",
        "    # split into vv and vh\n",
        "    df_vv = df_stats[df_stats['layer'].str.contains('VV')].reset_index().rename(columns={'layer': 'VV-pol'})\n",
        "    df_vh = df_stats[df_stats['layer'].str.contains('VH')].reset_index().rename(columns={'layer': 'VH-pol'})\n",
        "    \n",
        "    # concat vv and vh columns\n",
        "    df_con = pd.concat([df_vv, df_vh], axis=1, ignore_index=True)\n",
        "    \n",
        "    # rename columns\n",
        "    df_con.columns = ['i_2', 'lc_class', 'VV', 'Pixel count', 'VV mean', 'VV SD', 'VV slope', 'VV amp', 'lc class code',\n",
        "                      'i_3', 'lc_class_2', 'VH', 'VH count', 'VH mean', 'VH SD', 'VH slope', 'VH amp', 'lc class code_2']\n",
        "    \n",
        "    # subset columns\n",
        "    df_con = df_con[['lc_class', 'lc class code', 'VV', 'Pixel count', 'VV mean', 'VV SD', 'VV slope', 'VV amp', \n",
        "                     'VH', 'VH mean', 'VH SD', 'VH slope', 'VH amp']]\n",
        "    \n",
        "     # rename the layer names \n",
        "    df_con['VV'] = df_con['VV'].str.replace('VV_gamma0flat', '{}'.format(key))\n",
        "    df_con['VH'] = df_con['VH'].str.replace('VH_gamma0flat', '{}'.format(key))\n",
        "    \n",
        "    df_con['VV'] = df_con['VV'].str.replace('VV_gamma0', 'Original')\n",
        "    df_con['VH'] = df_con['VH'].str.replace('VH_gamma0', 'Original')\n",
        "    \n",
        "    # merge to existent dfs\n",
        "    df_con_merged.append(df_con)\n",
        "\n",
        "\n",
        "# bring all the data together\n",
        "appended_data = pd.concat(df_con_merged)\n",
        "\n",
        "# exclude marginal classes\n",
        "appended_data = appended_data[appended_data['Pixel count'] >= 100000]\n",
        "appended_data = appended_data[['lc_class', 'VV', 'VV mean', 'VV SD', 'VV slope', 'VV amp', 'VH mean', 'VH SD', 'VH slope', 'VH amp']]\n",
        "    \n",
        "\n",
        "# rename model names\n",
        "appended_data['VV'] = appended_data['VV'].str.replace('Bluesky_DTM_Clip_volume_buf_0', 'Model I')\n",
        "appended_data['VV'] = appended_data['VV'].str.replace('Bluesky_DTM_Clip_surface_buf_0', 'Model II')\n",
        "\n",
        "# remove double entries for original data and reindex \n",
        "appended_data.drop_duplicates(subset=['lc_class', 'VV'], keep='first', inplace=True)\n",
        "appended_data = appended_data.set_index(['lc_class', 'VV']).sort_values(['lc_class', 'VV'], ascending=False)\n",
        "\n",
        "# rename model column\n",
        "appended_data"
      ],
      "metadata": {
        "id": "wfLZEE01tR7B",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 318
        },
        "outputId": "976b5431-8d28-416b-f572-83241f0b8206"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "                                  VV mean     VV SD  VV slope    VV amp  \\\n",
              "lc_class              VV                                                  \n",
              "Improved Grassland    Original -10.463581  2.532052  0.094443  0.312362   \n",
              "                      Model II -10.473579  2.517432  0.035094  0.205654   \n",
              "                      Model I  -10.447763  2.518716  0.002001  0.203516   \n",
              "Broad-leaved forest   Original  -7.802244  2.024761  0.063566  0.476232   \n",
              "                      Model II  -7.918717  2.005812 -0.011170  0.318725   \n",
              "                      Model I   -7.868077  2.025835 -0.041572  0.394604   \n",
              "Arable & Horticulture Original  -7.860619  2.799963  0.079114  0.249216   \n",
              "                      Model II  -7.859532  2.796416  0.016141  0.208243   \n",
              "                      Model I   -7.846697  2.796992 -0.023514  0.213788   \n",
              "\n",
              "                                  VH mean     VH SD  VH slope    VH amp  \n",
              "lc_class              VV                                                 \n",
              "Improved Grassland    Original -17.335150  2.747841  0.094500  0.334236  \n",
              "                      Model II -17.345149  2.732126  0.035151  0.209755  \n",
              "                      Model I  -17.319332  2.733717  0.002058  0.190195  \n",
              "Broad-leaved forest   Original -13.106605  2.177710  0.064321  0.504687  \n",
              "                      Model II -13.223081  2.155194 -0.010410  0.337086  \n",
              "                      Model I  -13.172441  2.175344 -0.040812  0.402204  \n",
              "Arable & Horticulture Original -15.994920  2.913152  0.138289  0.356981  \n",
              "                      Model II -15.993827  2.906900  0.075229  0.279769  \n",
              "                      Model I  -15.980993  2.906576  0.035574  0.249910  "
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-fec10442-71f6-4419-a0bd-eeb3f8ad5207\">\n",
              "    <div class=\"colab-df-container\">\n",
              "      <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th>VV mean</th>\n",
              "      <th>VV SD</th>\n",
              "      <th>VV slope</th>\n",
              "      <th>VV amp</th>\n",
              "      <th>VH mean</th>\n",
              "      <th>VH SD</th>\n",
              "      <th>VH slope</th>\n",
              "      <th>VH amp</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>lc_class</th>\n",
              "      <th>VV</th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th rowspan=\"3\" valign=\"top\">Improved Grassland</th>\n",
              "      <th>Original</th>\n",
              "      <td>-10.463581</td>\n",
              "      <td>2.532052</td>\n",
              "      <td>0.094443</td>\n",
              "      <td>0.312362</td>\n",
              "      <td>-17.335150</td>\n",
              "      <td>2.747841</td>\n",
              "      <td>0.094500</td>\n",
              "      <td>0.334236</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Model II</th>\n",
              "      <td>-10.473579</td>\n",
              "      <td>2.517432</td>\n",
              "      <td>0.035094</td>\n",
              "      <td>0.205654</td>\n",
              "      <td>-17.345149</td>\n",
              "      <td>2.732126</td>\n",
              "      <td>0.035151</td>\n",
              "      <td>0.209755</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Model I</th>\n",
              "      <td>-10.447763</td>\n",
              "      <td>2.518716</td>\n",
              "      <td>0.002001</td>\n",
              "      <td>0.203516</td>\n",
              "      <td>-17.319332</td>\n",
              "      <td>2.733717</td>\n",
              "      <td>0.002058</td>\n",
              "      <td>0.190195</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th rowspan=\"3\" valign=\"top\">Broad-leaved forest</th>\n",
              "      <th>Original</th>\n",
              "      <td>-7.802244</td>\n",
              "      <td>2.024761</td>\n",
              "      <td>0.063566</td>\n",
              "      <td>0.476232</td>\n",
              "      <td>-13.106605</td>\n",
              "      <td>2.177710</td>\n",
              "      <td>0.064321</td>\n",
              "      <td>0.504687</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Model II</th>\n",
              "      <td>-7.918717</td>\n",
              "      <td>2.005812</td>\n",
              "      <td>-0.011170</td>\n",
              "      <td>0.318725</td>\n",
              "      <td>-13.223081</td>\n",
              "      <td>2.155194</td>\n",
              "      <td>-0.010410</td>\n",
              "      <td>0.337086</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Model I</th>\n",
              "      <td>-7.868077</td>\n",
              "      <td>2.025835</td>\n",
              "      <td>-0.041572</td>\n",
              "      <td>0.394604</td>\n",
              "      <td>-13.172441</td>\n",
              "      <td>2.175344</td>\n",
              "      <td>-0.040812</td>\n",
              "      <td>0.402204</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th rowspan=\"3\" valign=\"top\">Arable &amp; Horticulture</th>\n",
              "      <th>Original</th>\n",
              "      <td>-7.860619</td>\n",
              "      <td>2.799963</td>\n",
              "      <td>0.079114</td>\n",
              "      <td>0.249216</td>\n",
              "      <td>-15.994920</td>\n",
              "      <td>2.913152</td>\n",
              "      <td>0.138289</td>\n",
              "      <td>0.356981</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Model II</th>\n",
              "      <td>-7.859532</td>\n",
              "      <td>2.796416</td>\n",
              "      <td>0.016141</td>\n",
              "      <td>0.208243</td>\n",
              "      <td>-15.993827</td>\n",
              "      <td>2.906900</td>\n",
              "      <td>0.075229</td>\n",
              "      <td>0.279769</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Model I</th>\n",
              "      <td>-7.846697</td>\n",
              "      <td>2.796992</td>\n",
              "      <td>-0.023514</td>\n",
              "      <td>0.213788</td>\n",
              "      <td>-15.980993</td>\n",
              "      <td>2.906576</td>\n",
              "      <td>0.035574</td>\n",
              "      <td>0.249910</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-fec10442-71f6-4419-a0bd-eeb3f8ad5207')\"\n",
              "              title=\"Convert this dataframe to an interactive table.\"\n",
              "              style=\"display:none;\">\n",
              "        \n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
              "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
              "  </svg>\n",
              "      </button>\n",
              "      \n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      flex-wrap:wrap;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "      <script>\n",
              "        const buttonEl =\n",
              "          document.querySelector('#df-fec10442-71f6-4419-a0bd-eeb3f8ad5207 button.colab-df-convert');\n",
              "        buttonEl.style.display =\n",
              "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "        async function convertToInteractive(key) {\n",
              "          const element = document.querySelector('#df-fec10442-71f6-4419-a0bd-eeb3f8ad5207');\n",
              "          const dataTable =\n",
              "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                     [key], {});\n",
              "          if (!dataTable) return;\n",
              "\n",
              "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "            + ' to learn more about interactive tables.';\n",
              "          element.innerHTML = '';\n",
              "          dataTable['output_type'] = 'display_data';\n",
              "          await google.colab.output.renderOutput(dataTable, element);\n",
              "          const docLink = document.createElement('div');\n",
              "          docLink.innerHTML = docLinkHtml;\n",
              "          element.appendChild(docLink);\n",
              "        }\n",
              "      </script>\n",
              "    </div>\n",
              "  </div>\n",
              "  "
            ]
          },
          "metadata": {},
          "execution_count": 11
        }
      ]
    }
  ]
}