Copie_de_radar_based_flood_mapping.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/galython13/CRCPython/blob/master/Copie_de_radar_based_flood_mapping.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oHPuuAfTW3_I"
      },
      "source": [
        "<img src=\"https://github.com/vhertel/radar-based-flood-mapping/blob/main/resources/header.png?raw=1\" width=\"1000\"/>\n",
        "\n",
        "\n",
        "# Radar-based Flood Mapping\n",
        "\n",
        "<img src=\"https://github.com/vhertel/radar-based-flood-mapping/blob/main/resources/example.png?raw=1\" width=\"1000\"/>\n",
        "\n",
        "***\n",
        "\n",
        "The objective of this [Recommended Practice](https://un-spider.org/advisory-support/recommended-practices) is to determine the extent of flooded areas. The usage of Synthetic Aperture Radar (SAR) satellite imagery for flood extent mapping constitutes a viable solution with fast image processing, providing near real-time flood information to relief agencies for supporting humanitarian action. The high data reliability as well as the absence of geographical constraints, such as site accessibility, emphasize the technology’s potential in the field.\n",
        "\n",
        "This Jupyter Notebook is optimized for the use with Google Colab. As a cloud computing-based environment, it takes advantage of external technical resources and thus allows this tool to be applied using devices with limited computing power, including phones and tablets, and in areas with scarce bandwidth.  \n",
        "The notebook covers the full processing chain from data query and download up to the export of a final flood mask product by utilizing open access Sentinel-1 SAR data. The tool's workflow follows the UN-SPIDER Recommended Practice on [Radar-based Flood Mapping](https://un-spider.org/advisory-support/recommended-practices/recommended-practice-radar-based-flood-mapping) and is illustrated in the chart below. After entering user specifications, Sentinel-1  data can directly be downloaded from the <a href=\"https://scihub.copernicus.eu/\">Copernicus Open Access Hub</a>. Subsequently, the data is processed and stored in a variety of output formats.\n",
        "\n",
        "<img src=\"https://github.com/vhertel/radar-based-flood-mapping/blob/main/resources/charts/chart0.png?raw=1\" width=\"1000\"/>\n",
        "\n",
        "***\n",
        "\n",
        "***File Structure***  \n",
        "The Jupyter Notebook creates a folder called *'radar-based-flood-mapping'* in the Google Drive. Sentinel-1 images need to be stored in a subfolder called *'input'*. If no image is provided, the subfolder will automatically be created when accessing and downloading data from the <a href=\"https://scihub.copernicus.eu/\">Copernicus Open Access Hub</a> through this tool. If an area of interest (AOI) file is available (supported formats: GeoJSON, SHP, KML, KMZ), it needs to be placed in a subfolder called *'AOI'*. If none are available, a dialogue will allow to manually draw and/or upload locally stored AOI files. For reasons of automatic file selection, it is recommended to place only one AOI file in the respective folder. However, if multiple files exist, GeoJSON files are prioritized followed by SHP, KML, and KMZ. The processed data is stored in a subfolder called *'output'*.  \n",
        "In order to run the tool with no user interaction, all inputs must be clearly defined. This means that the *'input'* subfolder must include one single Sentinel-1 image and the *'AOI'* subfolder one single AOI file. All other scenarios do require manual interaction such as downloading data or defining an AOI.\n",
        "\n",
        "***Limitations***  \n",
        "Difficulties in detecting flooded vegetation and floods in urban areas due to double bounce backscatter. If water and non-water are very unequally distributed in the image, the histogram might not have a clear local minimum, leading to incorrect results in the automatic binarization process.\n",
        "\n",
        "***\n",
        "\n",
        "## Initialization\n",
        "\n",
        "The Jupyter Notebook takes advantage of the <a href=\"https://step.esa.int/docs/v6.0/apidoc/engine/\">ESA SNAP API Engine</a> using the SNAP-Python interface <i>snappy</i>. The installation and configuration procedure is included in the initialization step of this tool which consequently might take a few minutes during the initial run of this Jupyter Notebook."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "id": "EuHpDPr5eSnw",
        "outputId": "9038ccb2-39c1-4cdb-8c41-354b10f3c14e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ],
      "source": [
        "#@title <font color=#1B7192> Click to run </font>  { display-mode: \"form\" }\n",
        "\n",
        "#####################################################\n",
        "################### CONFIGURATION ###################\n",
        "#####################################################\n",
        "\n",
        "# mounnt Google Drive\n",
        "import os                                     # data access\n",
        "import google.colab                           # Google Colab\n",
        "import time                                   # time assessment\n",
        "import sys\n",
        "if not os.path.isdir('/content/drive'):\n",
        "    google.colab.drive.mount('/content/drive')\n",
        "\n",
        "try:\n",
        "    import snappy                             # SNAP Python interface\n",
        "    import jpy                                # Python-Java bridge\n",
        "except:\n",
        "    with google.colab.output.use_tags('snappy'):\n",
        "        sys.stdout.write('\\nPreparing conda environment...\\n')\n",
        "        sys.stdout.flush()\n",
        "        !pip install -q condacolab &> /dev/null\n",
        "        import condacolab\n",
        "        condacolab.install_miniconda()\n",
        "        sys.stdout.write('\\nInstalling snappy package...\\n')\n",
        "        sys.stdout.flush()\n",
        "        !conda install -c terradue -c conda-forge snap=8.0.0 &> /dev/null\n",
        "    google.colab.output.clear(output_tags='snappy')\n",
        "try:\n",
        "    import geopandas                          # data analysis and manipulation\n",
        "except:\n",
        "    with google.colab.output.use_tags('geopandas'):\n",
        "        sys.stdout.write('\\nInstalling geopandas package...\\n')\n",
        "        sys.stdout.flush()\n",
        "        !pip install geopandas &> /dev/null\n",
        "    google.colab.output.clear(output_tags='geopandas')\n",
        "try:\n",
        "    from sentinelsat.sentinel import SentinelAPI, read_geojson, geojson_to_wkt  # interface to Open Access Hub\n",
        "except:\n",
        "    with google.colab.output.use_tags('sentinelsat'):\n",
        "        sys.stdout.write('\\nInstalling sentinelsat package...\\n')\n",
        "        sys.stdout.flush()\n",
        "        !pip install sentinelsat &> /dev/null\n",
        "    google.colab.output.clear(output_tags='sentinelsat')\n",
        "try:\n",
        "    from ipyfilechooser import FileChooser    # file chooser widget\n",
        "except:\n",
        "    with google.colab.output.use_tags('ipyfilechooser'):\n",
        "        sys.stdout.write('\\nInstalling ipyfilechooser package...\\n')\n",
        "        sys.stdout.flush()\n",
        "        !pip install ipyfilechooser &> /dev/null\n",
        "    google.colab.output.clear(output_tags='ipyfilechooser')\n",
        "import ipywidgets\n",
        "if int(ipywidgets.__version__.split('.')[0]) >= 8:\n",
        "    !pip install \"ipywidgets>=7,<8\"  &> /dev/null\n",
        "\n",
        "# status update\n",
        "with google.colab.output.use_tags('initialization'):\n",
        "    sys.stdout.write('\\nInitialization successful!')\n",
        "    sys.stdout.flush()\n",
        "time.sleep(1)\n",
        "google.colab.output.clear(output_tags='initialization')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h0y2vRcGeJVB"
      },
      "source": [
        "## User Input\n",
        "\n",
        "<img src=\"https://github.com/vhertel/radar-based-flood-mapping/blob/main/resources/charts/chart1.png?raw=1\" width=\"1000\"/>\n",
        "\n",
        "Please specify in the code cell below i) the polarisation to be processed, ii) whether data shall be downloaded from the <a href=\"https://scihub.copernicus.eu/\">Copernicus Open Access Hub</a>, and iii) whether intermediate results should be plotted during the process. This section also loads relevant Python modules for the following analysis and initializes basic functionalities."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oX0CNEX7W3_P"
      },
      "outputs": [],
      "source": [
        "#@title <font color=#1B7192> Click to run </font>  { display-mode: \"form\" }\n",
        "\n",
        "####################################################\n",
        "#################### USER INPUT ####################\n",
        "####################################################\n",
        "\n",
        "# polarisations to be processed\n",
        "Polarisation  = 'VH'                    #@param [\"VH\", \"VV\", \"both\"]\n",
        "\n",
        "DownloadImage = True                   #@param {type:\"boolean\"}\n",
        "\n",
        "# show intermediate results if set to 'True'\n",
        "PlotResults   = True                    #@param {type:\"boolean\"}\n",
        "\n",
        "#####################################################\n",
        "###################### IMPORTS ######################\n",
        "#####################################################\n",
        "\n",
        "# MODULE                                      # DESCRIPTION\n",
        "import sys\n",
        "try:\n",
        "    import snappy                         # SNAP Python interface\n",
        "    import jpy                            # Python-Java bridge\n",
        "    import geopandas                          # data analysis and manipulation\n",
        "    from sentinelsat.sentinel import SentinelAPI, read_geojson, geojson_to_wkt  # interface to Open Access Hub\n",
        "    from ipyfilechooser import FileChooser    # file chooser widget\n",
        "except:\n",
        "    sys.exit('\\nPlease run initialization cell above first.')\n",
        "import os                                     # data access\n",
        "import google.colab                           # Google Colab\n",
        "import matplotlib.pyplot as plt               # create visualizations\n",
        "import numpy as np                            # scientific comupting\n",
        "import json                                   # JSON encoder and decoder\n",
        "import glob                                   # data access\n",
        "import ipywidgets                             # interactive UI controls\n",
        "import time                                   # time assessment\n",
        "import shutil                                 # file operations\n",
        "import folium                                 # visualization\n",
        "from folium import plugins                    # visualization\n",
        "from folium.plugins import MiniMap, Draw, Search # visualization\n",
        "import skimage.filters                        # threshold calculation\n",
        "import functools                              # higher-order functions and operations\n",
        "from datetime import date                     # dates, times and intervalls\n",
        "from IPython.display import display           # visualization\n",
        "from osgeo import ogr, gdal, osr              # data conversion\n",
        "from zipfile import ZipFile                   # file management\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "####################################################\n",
        "############### FUNCTION DEFINITIONS ###############\n",
        "####################################################\n",
        "\n",
        "def getAOI(path):\n",
        "    try:\n",
        "        file = readJSONFromAOI(path)\n",
        "    except:\n",
        "        print('No AOI file was found. Please draw and download the area of interest by clicking the Export-button')\n",
        "        print('inside the map or directly upload a locally stored AOI file using the dialogue below the map.\\n')\n",
        "        # create map\n",
        "        f = folium.Figure(height=500)\n",
        "        m = folium.Map(location=[0, 0], zoom_start=2.5, control_scale=True).add_to(f)\n",
        "        # add custom basemap\n",
        "        basemaps['Google Satellite Hybrid'].add_to(m)\n",
        "        # add a layer control panel to the map\n",
        "        m.add_child(folium.LayerControl())\n",
        "        # add minimap\n",
        "        m.add_child(MiniMap(tile_layer=basemaps['Google Satellite'], position='bottomright'))\n",
        "        # add draw control\n",
        "        draw = Draw(export=True, filename='AOI_manual_%s.geojson' % str(date.today()), draw_options={'polyline': False, 'circle': False, 'marker': False, 'circlemarker': False})\n",
        "        draw.add_to(m)\n",
        "        # display map\n",
        "        updater = display(f, display_id='m')\n",
        "        print('\\n')\n",
        "        # open upload section\n",
        "        os.chdir('/content')\n",
        "        uploaded = google.colab.files.upload()\n",
        "        for fn in uploaded.keys():\n",
        "            # copy uploaded file to GDrive folder\n",
        "            aoi_path = os.path.join(directory, 'AOI')\n",
        "            if not os.path.isdir(aoi_path):\n",
        "                os.mkdir(aoi_path)\n",
        "            shutil.copy2('/content/%s' % fn, aoi_path)\n",
        "            # remove original file\n",
        "            os.remove('/content/%s' % fn)\n",
        "            file_path = '%s/%s' % (aoi_path, fn)\n",
        "        file = readJSONFromAOI(aoi_path)\n",
        "\n",
        "    return file\n",
        "\n",
        "\n",
        "\n",
        "# Function looks for AOI file, converts to GeoJSON if not given and returns path to GeoJSON\n",
        "def readJSONFromAOI(path):\n",
        "    # check for GeoJSON file in 'AOI' subfolder\n",
        "    if len(glob.glob('%s/*.geojson' % path)) == 1:\n",
        "        file = glob.glob('%s/*.geojson' % path)[0]\n",
        "    elif len(glob.glob('%s/*.json' % path)) == 1:\n",
        "        file = glob.glob('%s/*.json' % path)[0]\n",
        "\n",
        "    # convert SHP to GeoJSON if no JSON is given\n",
        "    elif len(glob.glob('%s/*.shp' % path)) == 1:\n",
        "        file_name = os.path.splitext(glob.glob('%s/*.shp' % path)[0])[0].split('/')[-1]\n",
        "        shp_file = geopandas.read_file(glob.glob('%s/*.shp' % path)[0])\n",
        "        shp_file.to_file('%s/%s.json' % (path, file_name), driver='GeoJSON')\n",
        "        file = glob.glob('%s/*.json' % path)[0]\n",
        "\n",
        "    # convert KML to GeoJSON if no JSON or SHP is given\n",
        "    elif len(glob.glob('%s/*.kml' % path)) == 1:\n",
        "        file_name = os.path.splitext(glob.glob('%s/*.kml' % path)[0])[0].split('/')[-1]\n",
        "        kml_file = gdal.OpenEx(glob.glob('%s/*.kml' % path)[0])\n",
        "        ds = gdal.VectorTranslate('%s/%s.json' % (path, file_name), kml_file, format='GeoJSON')\n",
        "        del ds\n",
        "        file = glob.glob('%s/*.json' % path)[0]\n",
        "\n",
        "    # convert KMZ to JSON if no JSON, SHP, or KML is given\n",
        "    elif len(glob.glob('%s/*.kmz' % path)) == 1:\n",
        "        # open KMZ file and extract data\n",
        "        with ZipFile(glob.glob('%s/*.kmz' % path)[0], 'r') as kmz:\n",
        "            folder = os.path.splitext(glob.glob('%s/*.kmz' % path)[0])[0]\n",
        "            kmz.extractall(folder)\n",
        "        # convert KML to GeoJSON if extracted folder contains one KML file\n",
        "        if len(glob.glob('%s/*.kml' % folder)) == 1:\n",
        "            kml_file = gdal.OpenEx(glob.glob('%s/*.kml' % folder)[0])\n",
        "            ds = gdal.VectorTranslate('%s/%s.json' % (path, folder.split('/')[-1]), kml_file, format='GeoJSON')\n",
        "            del ds\n",
        "            file = glob.glob('%s/*.json' % path)[0]\n",
        "            # remove unzipped KMZ directory and data\n",
        "            shutil.rmtree(folder)\n",
        "    # allow to upload AOI file or manually draw AOI if no file was found\n",
        "    else:\n",
        "        raise FileNotFoundError\n",
        "\n",
        "    return file\n",
        "\n",
        "\n",
        "# plot band and histogram of 'Band'-type input and threshold\n",
        "# SNAP API: https://step.esa.int/docs/v6.0/apidoc/engine/\n",
        "def plotBand(band, threshold, binary=False):\n",
        "    # color stretch\n",
        "    vmin, vmax = 0, 1\n",
        "    # read pixel values\n",
        "    w = band.getRasterWidth()\n",
        "    h = band.getRasterHeight()\n",
        "    band_data = np.zeros(w * h, np.float32)\n",
        "    band.readPixels(0, 0, w, h, band_data)\n",
        "    band_data.shape = h, w\n",
        "    # color stretch\n",
        "    if binary:\n",
        "        cmap = plt.get_cmap('binary')\n",
        "    else:\n",
        "        vmin = np.percentile(band_data, 2.5)\n",
        "        vmax = np.percentile(band_data, 97.5)\n",
        "        cmap = plt.get_cmap('gray')\n",
        "    # plot band\n",
        "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
        "    ax1.imshow(band_data, cmap=cmap, vmin=vmin, vmax=vmax)\n",
        "    ax1.set_title(band.getName())\n",
        "    # plot histogram\n",
        "    band_data.shape = h * w\n",
        "    ax2.hist(np.asarray(band_data[band_data != 0], dtype='float'), bins=2048)\n",
        "    ax2.axvline(x=threshold, color='r')\n",
        "    ax2.set_title('Histogram: %s' % band.getName())\n",
        "\n",
        "    for ax in fig.get_axes():\n",
        "        ax.label_outer()\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "####################################################\n",
        "####################### CODE #######################\n",
        "####################################################\n",
        "\n",
        "# set working directory\n",
        "directory = '/content/drive/MyDrive/radar-based-flood-mapping'\n",
        "if not os.path.isdir(directory):\n",
        "    os.mkdir(directory)\n",
        "\n",
        "# Add custom basemap to folium\n",
        "basemaps = {\n",
        "    'Google Satellite Hybrid': folium.TileLayer(\n",
        "        tiles = 'https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}',\n",
        "        attr = 'Google',\n",
        "        name = 'Google Satellite',\n",
        "        overlay = True,\n",
        "        control = True,\n",
        "        show = False\n",
        "    ),\n",
        "    'Google Satellite': folium.TileLayer(\n",
        "        tiles = 'https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',\n",
        "        attr = 'Google',\n",
        "        name = 'Google Satellite',\n",
        "        overlay = True,\n",
        "        control = True,\n",
        "        show = False\n",
        "    )\n",
        "}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7aSEaXF1W3_S"
      },
      "source": [
        "## Download Image\n",
        "\n",
        "<img src=\"https://github.com/vhertel/radar-based-flood-mapping/blob/main/resources/charts/chart2.png?raw=1\" width=\"1000\"/>\n",
        "\n",
        "This section allows interactive data access and download from the <a href=\"https://scihub.copernicus.eu/\">Copernicus Open Access Hub</a>. It requires respective login details and the desired sensing period. If an AOI file is given in the *'AOI'* subfolder, the tool searches and displays available Sentinel-1 images accordingly. If no AOI file is provided, an interactive map allows to draw and download the area of interest by clicking the *'Export'*-button inside the map or to directly upload a locally stored AOI file. When hovering over a Sentinel-1 image, the tile index and ingestion dates are shown. The table below summarizes information on all available tiles and allows the download. The data is stored in the automatically created *'input'* subfolder. The Open Access Hub maintains an online archive of at least the latest year of products for immediate download. Access to previous products that are no longer available online will automatically trigger the retrieval from the long term archives. The actual download can be initiated by the user once the data are restored (within 24 hours)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wvsKF_MmW3_U"
      },
      "outputs": [],
      "source": [
        "#@title <font color=#1B7192> Click to run </font>  { display-mode: \"form\" }\n",
        "\n",
        "####################################################\n",
        "#################### USER INPUT ####################\n",
        "####################################################\n",
        "\n",
        "Username      = ''              #@param {type:\"string\"}\n",
        "Password      = ''              #@param {type:\"string\"}\n",
        "SensingPeriod_Start  = '2023-01-01'     #@param {type:\"date\"}\n",
        "SensingPeriod_Stop   = '2023-01-08'     #@param {type:\"date\"}\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "####################################################\n",
        "############### FUNCTION DEFINITIONS ###############\n",
        "####################################################\n",
        "\n",
        "# search for and display available Sentinel-1 tiles\n",
        "def queri(footprint):\n",
        "    # print status\n",
        "    with google.colab.output.use_tags('loading'):\n",
        "        sys.stdout.write('\\nLoading...')\n",
        "        sys.stdout.flush()\n",
        "    # search Copernicus Open Access Hub for products with regard to input footprint and sensing period\n",
        "    period = (date(int(SensingPeriod_Start.split('-')[0]), int(SensingPeriod_Start.split('-')[1]), int(SensingPeriod_Start.split('-')[2])),\n",
        "              date(int(SensingPeriod_Stop.split('-')[0]), int(SensingPeriod_Stop.split('-')[1]), int(SensingPeriod_Stop.split('-')[2])))\n",
        "    try:\n",
        "        products = api.query(footprint, date=period, platformname='Sentinel-1', producttype='GRD')\n",
        "        print('Successfully connected to Copernicus Open Access Hub.\\n', flush=True)\n",
        "    except:\n",
        "        sys.exit('\\nLogin data not valid. Please change username and/or password.')\n",
        "    # convert to GeoJSON for plot\n",
        "    products_json = api.to_geojson(products)\n",
        "    # raise warning that no image is available in given sensing period\n",
        "    if not products_json['features']:\n",
        "        sys.exit('\\nNo Sentinel-1 images available. Please change sensing period in user input section.')\n",
        "    # convert to dataframe for table visualization\n",
        "    products_df = api.to_dataframe(products)\n",
        "    # adds index to dataframe\n",
        "    indices = []\n",
        "    for i in range (1, len(products_df.index)+1):\n",
        "        indices.append('Tile %d' % i)\n",
        "        products_json.features[i-1].properties['index'] = ' Tile %d' % i\n",
        "    products_df.insert(0, 'index', indices, True)\n",
        "\n",
        "    # load products in folium framework for visualization\n",
        "    s1_tiles = folium.GeoJson(\n",
        "        products_json,\n",
        "        name='S1 tiles',\n",
        "        show=True,\n",
        "        style_function=lambda feature: {'fillColor': 'royalblue', 'fillOpacity' : 0.2},\n",
        "        highlight_function=lambda x: {'fillOpacity' : 0.4},\n",
        "        tooltip=folium.features.GeoJsonTooltip(\n",
        "            fields=['index', 'beginposition'],\n",
        "            aliases=['Index:','Date:'],\n",
        "        ),\n",
        "    ).add_to(m)\n",
        "    # add custom basemap\n",
        "    basemaps['Google Satellite Hybrid'].add_to(m)\n",
        "    # add a layer control panel to the map\n",
        "    m.add_child(folium.LayerControl())\n",
        "    # clear status update\n",
        "    google.colab.output.clear(output_tags='loading')\n",
        "    # update map\n",
        "    updater.update(m)\n",
        "\n",
        "    # print table with download buttons\n",
        "    grid = ipywidgets.GridspecLayout(len(products_df.index)+1, 5)\n",
        "    grid[0,0] = ipywidgets.HTML('<h3>Index</h3>')\n",
        "    grid[0,1] = ipywidgets.HTML('<h3>Date</h3>')\n",
        "    grid[0,2] = ipywidgets.HTML('<h3>Polarisation</h3>')\n",
        "    grid[0,3] = ipywidgets.HTML('<h3>Size</h3>')\n",
        "    for i in range(len(products_df.index)):\n",
        "        grid[i+1,0] = ipywidgets.Label(products_df['index'][i])\n",
        "        grid[i+1,1] = ipywidgets.Label(str(products_df['beginposition'][i]))\n",
        "        grid[i+1,2] = ipywidgets.Label(products_df['polarisationmode'][i])\n",
        "        grid[i+1,3] = ipywidgets.Label(products_df['size'][i])\n",
        "        grid[i+1,4] = ipywidgets.Button(description = 'Download')\n",
        "        grid[i+1,4].on_click(functools.partial(on_downloadButton_clicked, tile_id=products_df.values[i][-1]))\n",
        "    display(grid)\n",
        "\n",
        "\n",
        "# download chosen Sentinel-1 tile in subfolder 'input'\n",
        "def on_downloadButton_clicked(b, tile_id):\n",
        "    # get product information\n",
        "    product_info = api.get_product_odata(tile_id)\n",
        "    # check whether product is available\n",
        "    if product_info['Online']:\n",
        "        # check if input folder exists, if not create input folder\n",
        "        input_path = os.path.join(directory, 'input')\n",
        "        if not os.path.isdir(input_path):\n",
        "            os.mkdir(input_path)\n",
        "        # change into 'input' subfolder for storing product\n",
        "        os.chdir(input_path)\n",
        "        # status update\n",
        "        print('\\nProduct %s is online. Starting download.' % tile_id, flush=True)\n",
        "        # download product\n",
        "        api.download(tile_id)\n",
        "        # change back to previous working directory\n",
        "        os.chdir(directory)\n",
        "    # error message when product is not available\n",
        "    else:\n",
        "        print('\\nProduct %s is not online. Must be requested manually.\\n' % tile_id, flush=True)\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "####################################################\n",
        "####################### CODE #######################\n",
        "####################################################\n",
        "\n",
        "# check user input whether image download is requested\n",
        "if DownloadImage:\n",
        "    # connect to the API\n",
        "    api = SentinelAPI(Username, Password, 'https://scihub.copernicus.eu/dhus')\n",
        "    # get path to AOI\n",
        "    file = getAOI('%s/AOI' % directory)\n",
        "    # open AOI GeoJSON file and store data\n",
        "    with open(file, 'r') as f:\n",
        "        data_json = json.load(f)\n",
        "    # define map center depending on GeoJSON internal structure\n",
        "    try:\n",
        "        # GeoJSON format if KMZ is given\n",
        "        center = [data_json['features'][0]['geometry']['coordinates'][0][0][0][1],\n",
        "                  data_json['features'][0]['geometry']['coordinates'][0][0][0][0]]\n",
        "    except:\n",
        "        # GeoJSON format if JSON or SHP is given\n",
        "        center = [data_json['features'][0]['geometry']['coordinates'][0][0][1],\n",
        "                  data_json['features'][0]['geometry']['coordinates'][0][0][0]]\n",
        "    # create map\n",
        "    f = folium.Figure(height=500)\n",
        "    m = folium.Map(location=center, zoom_start=6, control_scale=True).add_to(f)\n",
        "    # add AOI to map\n",
        "    folium.GeoJson(file, name='AOI', style_function = lambda x: {'color':'green'}).add_to(m)\n",
        "    footprint = geojson_to_wkt(data_json)\n",
        "    updater = display(f, display_id='m')\n",
        "\n",
        "    # search for available Sentinel-1 tiles\n",
        "    queri(footprint)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xoLWzaiXW3_V"
      },
      "source": [
        "## Processing\n",
        "\n",
        "<img src=\"https://github.com/vhertel/radar-based-flood-mapping/blob/main/resources/charts/chart3.png?raw=1\" width=\"1000\"/>\n",
        "\n",
        "If more than one Sentinel-1 image exists in the *'input'* subfolder, the user can select which one is to be used for the processing. The subset is generated according to the AOI file in the *'AOI'* subfolder. If no AOI file is provided, an interactive map allows to draw and download the area of interest by clicking the *'Export'*-button inside the map or to directly upload a locally stored AOI file. Subsequently, the following processing steps are performed:\n",
        "\n",
        "1. ***Apply Orbit File***: The orbit file provides accurate satellite position and velocity information. Based on this information, the orbit state vectors in the abstract metadata of the product are updated. The precise orbit files are available days-to-weeks after the generation of the product. Since this is an optional processing step, the tool will continue the workflow in case the orbit file is not yet available to allow rapid mapping applications.\n",
        "\n",
        "\n",
        "2. ***Thermal Noise Removal***: Thermal noise correction is applied to Sentinel-1 Level-1 GRD products which have not already been corrected.\n",
        "\n",
        "\n",
        "3. ***Radiometric Calibration***: The objective of SAR calibration is to provide imagery in which the pixel values can be directly related to the radar backscatter of the scene. Though uncalibrated SAR imagery is sufficient for qualitative use, calibrated SAR images are essential to the quantitative use of SAR data.\n",
        "\n",
        "\n",
        "4. ***Speckle Filtering***: SAR images have inherent texturing called speckles which degrade the quality of the image and make interpretation of features more difficult. Speckles are caused by random constructive and destructive interference of the de-phased but coherent return waves scattered by the elementary scatter within each resolution cell. Speckle noise reduction can be applied either by spatial filtering or multilook processing. A Lee filter with an X, Y size of 5, 5 is used in this step.\n",
        "\n",
        "\n",
        "5. ***Terrain Correction***: Due to topographical variations of a scene and the tilt of the satellite sensor, distances can be distorted in the SAR images. Data which is not directly directed towards the sensor’s Nadir location will have some distortion. Therefore, terrain corrections are intended to compensate for these distortions to allow a realistic geometric representation in the image.\n",
        "\n",
        "\n",
        "6. ***Binarization***: In order to obtain a binary flood mask, the histogram is analyzed to separate water from non-water pixels. Due to the side-looking geometry of SAR sensors and the comparably smooth surface of water, only a very small proportion of backscatter is reflected back to the sensor leading to comparably low pixel values in the histogram. The threshold used for separation is automatically calculated using <a href=\"https://scikit-image.org/\">scikit-image</a> implementations and a combined use of the <a href=\"https://doi.org/10.1111/j.1749-6632.1965.tb11715.x\">minimum method</a> and <a href=\"https://www.semanticscholar.org/paper/A-threshold-selection-method-from-gray-level-Otsu/1d4816c612e38dac86f2149af667a5581686cdef?p2df\">Otsu's method</a>. The <a href=\"http://due.esrin.esa.int/page_globcover.php\">GlobCover</a> layer of the European Space Agency is used to mask out permanent water bodies.\n",
        "\n",
        "\n",
        "7. ***Speckle Filtering***: A Median filter with an X, Y size of 7, 7 is used in this step."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "q3yEG-4jW3_V"
      },
      "outputs": [],
      "source": [
        "#@title <font color=#1B7192> Click to run </font>  { display-mode: \"form\" }\n",
        "\n",
        "####################################################\n",
        "############### FUNCTION DEFINITIONS ###############\n",
        "####################################################\n",
        "\n",
        "def applySubset(path):\n",
        "    # set correct path of input file and create S1 product\n",
        "    if len(files) == 1:\n",
        "        file_path = path\n",
        "    else:\n",
        "        file_path = path.selected\n",
        "    S1_source = snappy.ProductIO.readProduct(file_path)\n",
        "\n",
        "    # read geographic coordinates from Sentinel-1 image meta data\n",
        "    meta_data = S1_source.getMetadataRoot().getElement('Abstracted_Metadata')\n",
        "    # refines center of map according to Sentinel-1 image\n",
        "    S1_center = (meta_data.getAttributeDouble('centre_lat'), meta_data.getAttributeDouble('centre_lon'))\n",
        "    # defines polygon illustrating Sentinel-1 image\n",
        "    polygon_geom = {\n",
        "      \"type\": \"FeatureCollection\",\n",
        "      \"features\":\n",
        "                [{\"type\": \"Feature\",\n",
        "                \"properties\": {},\n",
        "                \"geometry\": {\"type\": \"Polygon\", \"coordinates\": [[[meta_data.getAttributeDouble('first_near_long'), meta_data.getAttributeDouble('first_near_lat')],\n",
        "                                                                [meta_data.getAttributeDouble('last_near_long'), meta_data.getAttributeDouble('last_near_lat')],\n",
        "                                                                [meta_data.getAttributeDouble('last_far_long'), meta_data.getAttributeDouble('last_far_lat')],\n",
        "                                                                [meta_data.getAttributeDouble('first_far_long'), meta_data.getAttributeDouble('first_far_lat')],\n",
        "                                                                [meta_data.getAttributeDouble('first_near_long'), meta_data.getAttributeDouble('first_near_lat')]]]}}]}\n",
        "\n",
        "    # get path to AOI\n",
        "    file = getAOI('%s/AOI' % directory)\n",
        "    # open GeoJSON file and store data\n",
        "    with open(file, 'r') as f:\n",
        "        data_json = json.load(f)\n",
        "    footprint = geojson_to_wkt(data_json)\n",
        "\n",
        "    # creates map\n",
        "    f = folium.Figure(height=500)\n",
        "    m = folium.Map(location = S1_center, zoom_start = 7.5, control_scale=True).add_to(f)\n",
        "    # add S1 tile to map\n",
        "    folium.GeoJson(polygon_geom, name='Sentinel-1 tile').add_to(m)\n",
        "    # add AOI to map\n",
        "    folium.GeoJson(file, name='AOI', style_function = lambda x: {'color':'green'}).add_to(m)\n",
        "    # add custom basemap\n",
        "    basemaps['Google Satellite Hybrid'].add_to(m)\n",
        "    # add a layer control panel to the map\n",
        "    m.add_child(folium.LayerControl())\n",
        "    # display map\n",
        "    updater = display(f, display_id='m')\n",
        "\n",
        "    # Subset operator\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('copyMetadata', True)\n",
        "    geom = snappy.WKTReader().read(footprint)\n",
        "    parameters.put('geoRegion', geom)\n",
        "    parameters.put('sourceBands', sourceBands)\n",
        "    S1_crop = snappy.GPF.createProduct('Subset', parameters, S1_source)\n",
        "    # status update\n",
        "    print('\\nSubset successfully generated.\\n', flush=True)\n",
        "\n",
        "    # run processing process\n",
        "    processing(S1_crop)\n",
        "\n",
        "\n",
        "# calculate and return threshold of 'Band'-type input\n",
        "# SNAP API: https://step.esa.int/docs/v6.0/apidoc/engine/\n",
        "def getThreshold(S1_band):\n",
        "    # read band\n",
        "    w = S1_band.getRasterWidth()\n",
        "    h = S1_band.getRasterHeight()\n",
        "    band_data = np.zeros(w * h, np.float32)\n",
        "    S1_band.readPixels(0, 0, w, h, band_data)\n",
        "    band_data.shape = h * w\n",
        "    # calculate threshold using Otsu method\n",
        "    threshold_otsu = skimage.filters.threshold_otsu(band_data)\n",
        "    # calculate threshold using minimum method\n",
        "    threshold_minimum = skimage.filters.threshold_minimum(band_data)\n",
        "    # get number of pixels for both thresholds\n",
        "    numPixOtsu = len(band_data[abs(band_data - threshold_otsu) < 0.1])\n",
        "    numPixMinimum = len(band_data[abs(band_data - threshold_minimum) < 0.1])\n",
        "\n",
        "    # if number of pixels at minimum threshold is less than 0.1% of number of pixels at Otsu threshold\n",
        "    if abs(numPixMinimum/numPixOtsu) < 0.001:\n",
        "        # adjust band data according\n",
        "        if threshold_otsu < threshold_minimum:\n",
        "            band_data = band_data[band_data < threshold_minimum]\n",
        "            threshold_minimum = skimage.filters.threshold_minimum(band_data)\n",
        "        else:\n",
        "            band_data = band_data[band_data > threshold_minimum]\n",
        "            threshold_minimum = skimage.filters.threshold_minimum(band_data)\n",
        "        numPixMinimum = len(band_data[abs(band_data - threshold_minimum) < 0.1])\n",
        "    # check for final threshold\n",
        "    if abs(numPixMinimum/numPixOtsu) < 0.001:\n",
        "        threshold = threshold_otsu\n",
        "    else:\n",
        "        threshold = threshold_minimum\n",
        "\n",
        "    return threshold\n",
        "\n",
        "\n",
        "# calculate binary mask of 'Product'-type intput with respect expression in string array\n",
        "def binarization(S1_product, expressions):\n",
        "\n",
        "    BandDescriptor = jpy.get_type('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor')\n",
        "    targetBands = jpy.array('org.esa.snap.core.gpf.common.BandMathsOp$BandDescriptor', len(expressions))\n",
        "\n",
        "    # loop through bands\n",
        "    for i in range(len(expressions)):\n",
        "        targetBand = BandDescriptor()\n",
        "        targetBand.name = '%s' % S1_product.getBandNames()[i]\n",
        "        targetBand.type = 'float32'\n",
        "        targetBand.expression = expressions[i]\n",
        "        targetBands[i] = targetBand\n",
        "\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('targetBands', targetBands)\n",
        "    mask = snappy.GPF.createProduct('BandMaths', parameters, S1_product)\n",
        "\n",
        "    return mask\n",
        "\n",
        "\n",
        "# processing steps\n",
        "def processing(S1_crop):\n",
        "    # Apply-Orbit-File operator\n",
        "    print('1. Apply Orbit File:          ', end='', flush=True)\n",
        "    start_time = time.time()\n",
        "    parameters = snappy.HashMap()\n",
        "    # continue with calculation in case no orbit file is available yet\n",
        "    parameters.put('continueOnFail', True)\n",
        "    S1_Orb = snappy.GPF.createProduct('Apply-Orbit-File', parameters, S1_crop)\n",
        "    print('--- %.2f  seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "    # ThermalNoiseRemoval operator\n",
        "    print('2. Thermal Noise Removal:     ', end='', flush=True)\n",
        "    start_time = time.time()\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('removeThermalNoise', True)\n",
        "    S1_Thm = snappy.GPF.createProduct('ThermalNoiseRemoval', parameters, S1_Orb)\n",
        "    print('--- %.2f  seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "    # Calibration operator\n",
        "    print('3. Radiometric Calibration:   ', end='', flush=True)\n",
        "    start_time = time.time()\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('outputSigmaBand', True)\n",
        "    S1_Cal = snappy.GPF.createProduct('Calibration', parameters, S1_Thm)\n",
        "    print('--- %.2f  seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "    # Speckle-Filter operator\n",
        "    print('4. Speckle Filtering:         ', end='', flush=True)\n",
        "    start_time = time.time()\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('filter', 'Lee')\n",
        "    parameters.put('filterSizeX', 5)\n",
        "    parameters.put('filterSizeY', 5)\n",
        "    S1_Spk = snappy.GPF.createProduct('Speckle-Filter', parameters, S1_Cal)\n",
        "    print('--- %.2f  seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "    # Conversion from linear to db operator\n",
        "    S1_Spk_db = snappy.GPF.createProduct('LinearToFromdB', snappy.HashMap(), S1_Spk)\n",
        "\n",
        "    # Terrain-Correction operator\n",
        "    print('5. Terrain Correction:        ', end='', flush=True)\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('demName', 'SRTM 1Sec HGT')\n",
        "    parameters.put('demResamplingMethod', 'BILINEAR_INTERPOLATION')\n",
        "    parameters.put('imgResamplingMethod', 'NEAREST_NEIGHBOUR')\n",
        "    parameters.put('pixelSpacingInMeter', 10.0)\n",
        "    parameters.put('nodataValueAtSea', False)\n",
        "    parameters.put('saveSelectedSourceBand', True)\n",
        "    S1_TC = snappy.GPF.createProduct('Terrain-Correction', parameters, S1_Spk_db)\n",
        "    print('--- %.2f  seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "    # Binarization\n",
        "    print('6. Binarization:              ', end='', flush=True)\n",
        "    start_time = time.time()\n",
        "    # add GlobCover band\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('landCoverNames', 'GlobCover')\n",
        "    GlobCover = snappy.GPF.createProduct('AddLandCover', parameters, S1_TC)\n",
        "    # empty string array for binarization band maths expression(s)\n",
        "    expressions = ['' for i in range(S1_TC.getNumBands())]\n",
        "    # empty array for threshold(s)\n",
        "    thresholds = np.zeros(S1_TC.getNumBands())\n",
        "    # loop through bands\n",
        "    for i in range(S1_TC.getNumBands()):\n",
        "        # calculate threshold of band and store in float array\n",
        "        # use S1_Spk_db product for performance reasons. S1_TC causes 0-values\n",
        "        # which distort histogram and thus threshold result\n",
        "        thresholds[i] = getThreshold(S1_Spk_db.getBandAt(i))\n",
        "        # formulate expression according to threshold and store in string array\n",
        "        expressions[i] = 'if (%s < %s && land_cover_GlobCover != 210) then 1 else NaN' % (S1_TC.getBandNames()[i], thresholds[i])\n",
        "    # do binarization\n",
        "    S1_floodMask = binarization(GlobCover, expressions)\n",
        "    print('--- %.2f seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "    # Speckle-Filter operator\n",
        "    print('7. Speckle Filtering:         ', end='', flush=True)\n",
        "    start_time = time.time()\n",
        "    parameters = snappy.HashMap()\n",
        "    parameters.put('filter', 'Median')\n",
        "    parameters.put('filterSizeX', 5)\n",
        "    parameters.put('filterSizeY', 5)\n",
        "    # define flood mask as global for later access\n",
        "    global S1_floodMask_Spk\n",
        "    S1_floodMask_Spk = snappy.GPF.createProduct('Speckle-Filter', parameters, S1_floodMask)\n",
        "    print('--- %.2f  seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "    # output\n",
        "    if PlotResults:\n",
        "        print('8. Plot:                      ', end='', flush=True)\n",
        "        start_time = time.time()\n",
        "        for i in range(S1_TC.getNumBands()):\n",
        "            plotBand(S1_TC.getBandAt(i), thresholds[i])\n",
        "        print('--- %.2f seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "####################################################\n",
        "####################### CODE #######################\n",
        "####################################################\n",
        "\n",
        "# filter required polarisation(s) and set output file name accordingly\n",
        "if Polarisation == 'both':\n",
        "    sourceBands = 'Amplitude_VH,Intensity_VH,Amplitude_VV,Intensity_VV'\n",
        "    output_extensions   = 'processed_VHVV'\n",
        "elif Polarisation == 'VH':\n",
        "    sourceBands = 'Amplitude_VH,Intensity_VH'\n",
        "    output_extensions   = 'processed_VH'\n",
        "elif Polarisation == 'VV':\n",
        "    sourceBands = 'Amplitude_VV,Intensity_VV'\n",
        "    output_extensions   = 'processed_VV'\n",
        "\n",
        "# path of Sentinel-1 .zip input file\n",
        "input_path = os.path.join(directory, 'input')\n",
        "# empty string array to store Sentinel-1 files in 'input' subfolder\n",
        "files = []\n",
        "# add files to list\n",
        "for file in glob.glob1(input_path, '*.zip'):\n",
        "    files.append(file)\n",
        "# select input file and start processing if there is only one available Sentinel-1 file\n",
        "if len(files) == 1:\n",
        "    input_name = files[0]\n",
        "    print('Selected:  %s\\n' % input_name, flush=True)\n",
        "    # apply subset according to JSON data\n",
        "    applySubset('%s/%s' % (input_path, input_name))\n",
        "# open dialogue to select input file if more or less than one is available\n",
        "else:\n",
        "    print('More or less than one Sentinel-1 files have been found. Please select.', flush=True)\n",
        "    fc = FileChooser(input_path)\n",
        "    fc.filter_pattern = '*.zip'\n",
        "    fc.register_callback(applySubset)\n",
        "    display(fc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MavzlDEFW3_W"
      },
      "source": [
        "## Data Export\n",
        "\n",
        "<img src=\"https://github.com/vhertel/radar-based-flood-mapping/blob/main/resources/charts/chart4.png?raw=1\" width=\"1000\"/>\n",
        "\n",
        "The processed flood mask is exported as GeoTIFF, SHP, KML, and GeoJSON and stored in the *'output'* subfolder. An interactive map shows the flood mask."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "og4yIuW3W3_W"
      },
      "outputs": [],
      "source": [
        "#@title <font color=#1B7192> Click to run </font>  { display-mode: \"form\" }\n",
        "\n",
        "####################################################\n",
        "####################### CODE #######################\n",
        "####################################################\n",
        "\n",
        "print('Exporting...\\n', flush=True)\n",
        "# check if output folders exists, if not create folders\n",
        "output_path = os.path.join(directory, 'output')\n",
        "if not os.path.isdir(output_path):\n",
        "    os.mkdir(output_path)\n",
        "GeoTIFF_path = os.path.join(output_path, 'GeoTIFF')\n",
        "if not os.path.isdir(GeoTIFF_path):\n",
        "    os.mkdir(GeoTIFF_path)\n",
        "SHP_path = os.path.join(output_path, 'SHP')\n",
        "if not os.path.isdir(SHP_path):\n",
        "    os.mkdir(SHP_path)\n",
        "KML_path = os.path.join(output_path, 'KML')\n",
        "if not os.path.isdir(KML_path):\n",
        "    os.mkdir(KML_path)\n",
        "GeoJSON_path = os.path.join(output_path, 'GeoJSON')\n",
        "if not os.path.isdir(GeoJSON_path):\n",
        "    os.mkdir(GeoJSON_path)\n",
        "# get file name if file chooser was used\n",
        "if len(files) != 1: input_name = fc.selected_filename\n",
        "\n",
        "# write output file as GeoTIFF\n",
        "print('1. GeoTIFF:                   ', end='', flush=True)\n",
        "start_time = time.time()\n",
        "snappy.ProductIO.writeProduct(S1_floodMask_Spk, '%s/%s_%s' % (GeoTIFF_path, os.path.splitext(input_name)[0], output_extensions), 'GeoTIFF')\n",
        "print('--- %.2f seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "# convert GeoTIFF to SHP\n",
        "print('2. SHP:                       ', end='', flush=True)\n",
        "start_time = time.time()\n",
        "# allow GDAL to throw Python exceptions\n",
        "gdal.UseExceptions()\n",
        "open_image = gdal.Open('%s/%s_%s.tif' % (GeoTIFF_path, os.path.splitext(input_name)[0], output_extensions))\n",
        "srs = osr.SpatialReference()\n",
        "srs.ImportFromWkt(open_image.GetProjectionRef())\n",
        "shp_driver = ogr.GetDriverByName('ESRI Shapefile')\n",
        "# empty string array for bands in GeoTIFF\n",
        "output_shp = ['' for i in range(open_image.RasterCount)]\n",
        "if open_image.RasterCount == 1:\n",
        "    output_shp[0] = '%s/%s_processed_%s' % (SHP_path, os.path.splitext(input_name)[0], Polarisation)\n",
        "else:\n",
        "    VH_SHP_path = os.path.join(SHP_path, 'VH')\n",
        "    if not os.path.isdir(VH_SHP_path):\n",
        "        os.mkdir(VH_SHP_path)\n",
        "    VV_SHP_path = os.path.join(SHP_path, 'VV')\n",
        "    if not os.path.isdir(VV_SHP_path):\n",
        "        os.mkdir(VV_SHP_path)\n",
        "    output_shp[0] = '%s/%s_processed_VH' % (VH_SHP_path, os.path.splitext(input_name)[0])\n",
        "    output_shp[1] = '%s/%s_processed_VV' % (VV_SHP_path, os.path.splitext(input_name)[0])\n",
        "# loops through bands in GeoTIFF\n",
        "for i in range(open_image.RasterCount):\n",
        "    input_band = open_image.GetRasterBand(i+1)\n",
        "    output_shapefile = shp_driver.CreateDataSource(output_shp[i] + '.shp')\n",
        "    new_shapefile = output_shapefile.CreateLayer(output_shp[i], srs=srs)\n",
        "    new_shapefile.CreateField(ogr.FieldDefn('DN', ogr.OFTInteger))\n",
        "    gdal.Polygonize(input_band, input_band.GetMaskBand(), new_shapefile, 0, [], callback=None)\n",
        "    # filters attributes with values other than 1 (sould be NaN or respective value)\n",
        "    new_shapefile.SetAttributeFilter('DN != 1')\n",
        "    for feat in new_shapefile:\n",
        "        new_shapefile.DeleteFeature(feat.GetFID())\n",
        "    new_shapefile.SyncToDisk()\n",
        "print('--- %.2f seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "# convert SHP to KML\n",
        "print('3. KML:                       ', end='', flush=True)\n",
        "start_time = time.time()\n",
        "if open_image.RasterCount == 1:\n",
        "    shp_file = gdal.OpenEx('%s/%s_processed_%s.shp' % (SHP_path, os.path.splitext(input_name)[0], Polarisation))\n",
        "    ds = gdal.VectorTranslate('%s/%s_processed_%s.kml' % (KML_path, os.path.splitext(input_name)[0], Polarisation), shp_file, format='KML')\n",
        "    del ds\n",
        "else:\n",
        "    shp_file_VH = gdal.OpenEx('%s/%s_processed_VH.shp' % (VH_SHP_path, os.path.splitext(input_name)[0]))\n",
        "    ds_VH = gdal.VectorTranslate('%s/%s_processed_VH.kml' % (KML_path, os.path.splitext(input_name)[0]), shp_file_VH, format='KML')\n",
        "    del ds_VH\n",
        "    shp_file_VV = gdal.OpenEx('%s/%s_processed_VV.shp' % (VV_SHP_path, os.path.splitext(input_name)[0]))\n",
        "    ds_VV = gdal.VectorTranslate('%s/%s_processed_VV.kml' % (KML_path, os.path.splitext(input_name)[0]), shp_file_VV, format='KML')\n",
        "    del ds_VV\n",
        "print('--- %.2f seconds ---' % (time.time() - start_time), flush=True)\n",
        "\n",
        "# convert SHP to GeoJSON\n",
        "print('4. GeoJSON:                   ', end='', flush=True)\n",
        "start_time = time.time()\n",
        "if open_image.RasterCount == 1:\n",
        "    shp_file = geopandas.read_file('%s/%s_processed_%s.shp' % (SHP_path, os.path.splitext(input_name)[0], Polarisation))\n",
        "    shp_file.to_file('%s/%s_processed_%s.json' % (GeoJSON_path, os.path.splitext(input_name)[0], Polarisation), driver='GeoJSON')\n",
        "else:\n",
        "    shp_file_VH = geopandas.read_file('%s/%s_processed_VH.shp' % (VH_SHP_path, os.path.splitext(input_name)[0]))\n",
        "    shp_file_VH.to_file('%s/%s_processed_VH.json' % (GeoJSON_path, os.path.splitext(input_name)[0]), driver='GeoJSON')\n",
        "    shp_file_VV = geopandas.read_file('%s/%s_processed_VV.shp' % (VV_SHP_path, os.path.splitext(input_name)[0]))\n",
        "    shp_file_VV.to_file('%s/%s_processed_VV.json' % (GeoJSON_path, os.path.splitext(input_name)[0]), driver='GeoJSON')\n",
        "print('--- %.2f seconds ---\\n' % (time.time() - start_time), flush=True)\n",
        "print('Files successfuly stored under %s.\\n' % output_path, flush=True)\n",
        "\n",
        "# plot results\n",
        "meta_data = S1_floodMask_Spk.getMetadataRoot().getElement('Abstracted_Metadata')\n",
        "S1_center = (meta_data.getAttributeDouble('centre_lat'), meta_data.getAttributeDouble('centre_lon'))\n",
        "f = folium.Figure(height=500)\n",
        "results_map = folium.Map(location = S1_center, zoom_start = 9, control_scale = True).add_to(f)\n",
        "if open_image.RasterCount == 1:\n",
        "    file = '%s/%s_processed_%s.json' % (GeoJSON_path, os.path.splitext(input_name)[0], Polarisation)\n",
        "    folium.GeoJson(file, name='Flood Mask %s' % Polarisation, style_function = lambda x: {'color':'blue', 'opacity':'1', 'fillColor':'blue', 'fillOpacity':'1', 'weight':'0.8'}).add_to(results_map)\n",
        "else:\n",
        "    file_VV = '%s/%s_processed_VV.json' % (GeoJSON_path, os.path.splitext(input_name)[0])\n",
        "    folium.GeoJson(file_VV, name='Flood Mask VV', style_function = lambda x: {'color':'red', 'opacity':'1', 'fillColor':'red', 'fillOpacity':'1', 'weight':'0.8'}).add_to(results_map)\n",
        "    file_VH = '%s/%s_processed_VH.json' % (GeoJSON_path, os.path.splitext(input_name)[0])\n",
        "    folium.GeoJson(file_VH, name='Flood Mask VH', style_function = lambda x: {'color':'blue', 'opacity':'1', 'fillColor':'blue', 'fillOpacity':'1', 'weight':'0.8'}).add_to(results_map)\n",
        "# add custom basemap\n",
        "basemaps['Google Satellite Hybrid'].add_to(results_map)\n",
        "# add a layer control panel to the map\n",
        "results_map.add_child(folium.LayerControl(collapsed=False))\n",
        "display(f)"
      ]
    }
  ],
  "metadata": {
    "accelerator": "TPU",
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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": {
        "height": "721px",
        "left": "82px",
        "top": "111.133px",
        "width": "237px"
      },
      "toc_section_display": true,
      "toc_window_display": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}