radar-based-flood-mapping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "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 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 file constitutes the directory of origin. Additional data is contained in subfolders. 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, an interactive map will allow to manually draw the area of interest. 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",
    "***Important***\n",
    "The Jupyter Notebook takes advantage of the <a href=\"https://step.esa.int/docs/v6.0/apidoc/engine/\">ESA SNAP API Engine</a> and requires installation of the SNAP-Python interface <i>snappy</i>. Click <a href=\"https://senbox.atlassian.net/wiki/spaces/SNAP/pages/24051781/Using+SNAP+in+your+Python+programs\">here</a> for further information. Furthermore, the <a href=\"https://jupyter-contrib-nbextensions.readthedocs.io/en/latest/install.html\">Jupyter Notebook Extensions</a> *Codefolding*, *ExecuteTime* and *Table of Contents (2)* are used for the most convene performance.\n",
    "\n",
    "***\n",
    "\n",
    "## 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> with respective sensing period and login details, and iii) whether intermediate results should be plotted during the process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# polarisations to be processed\n",
    "polarisations = 'VH'                              # 'VH', 'VV', 'both'\n",
    "\n",
    "# download image from Copernicus Open Access Hub\n",
    "download = {\n",
    "    'imageDownload'     : True,                   # 'True', 'False'\n",
    "    'period_start'      : [2021, 1, 1],           # format: [Year, Month, Day]\n",
    "    'period_stop'       : [2021, 1, 8],           # format: [Year, Month, Day]\n",
    "    'username'          : 'username',             # username for login\n",
    "    'password'          : 'password'              # password for login\n",
    "}\n",
    "\n",
    "# show intermediate results if set to 'True'\n",
    "plotResoluts = True                               # 'True', 'False'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization\n",
    "\n",
    "This section loads relevant Python modules for the following analysis and initializes basic functionalities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     37,
     85
    ]
   },
   "outputs": [],
   "source": [
    "# Click to run\n",
    "\n",
    "#####################################################\n",
    "###################### IMPORTS ######################\n",
    "#####################################################\n",
    "\n",
    "# MODULE                                      # DESCRIPTION\n",
    "import sys\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 os                                     # data access\n",
    "import ipywidgets                             # interactive UI controls\n",
    "import time                                   # time assessment\n",
    "import shutil                                 # file operations\n",
    "import ipyleaflet                             # visualization\n",
    "import geopandas                              # data analysis and manipulation\n",
    "import snappy                                 # SNAP Python interface\n",
    "import jpy                                    # Python-Java bridge\n",
    "import skimage.filters                        # threshold calculation\n",
    "import functools                              # higher-order functions and operations\n",
    "from ipyfilechooser import FileChooser        # file chooser widget\n",
    "from sentinelsat.sentinel import SentinelAPI, read_geojson, geojson_to_wkt  # interface to Open Access Hub\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",
    "# Function looks for AOI file, converts to GeoJSON if not given and returns 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",
    "    # open JSON file and store data\n",
    "    with open(file, 'r') as f:\n",
    "        data_json = json.load(f)\n",
    "\n",
    "    return data_json\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",
    "# get current working directory\n",
    "directory = os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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>. 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, the search bar on the left side of the interactive map can be used to find the desired region. The AOI can be then be selected and manipulated manually by using the drawing tool. Clicking the Search button below the map will load available images. If multiple AOIs are drawn, only the last one is considered. When hovering over a Sentinel-1 image, the tile index and ingestion date 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": {
    "code_folding": [
     0,
     7,
     57,
     65,
     79
    ]
   },
   "outputs": [],
   "source": [
    "# Click to run\n",
    "\n",
    "####################################################\n",
    "############### FUNCTION DEFINITIONS ###############\n",
    "####################################################\n",
    "\n",
    "# search for and display available Sentinel-1 tiles\n",
    "def queri(footprint):\n",
    "    # print status\n",
    "    print('Loading...', flush=True)\n",
    "    # search Copernicus Open Access Hub for products with regard to input footprint and sensing period\n",
    "    try:\n",
    "        products = api.query(footprint,\n",
    "                             date = (date(download['period_start'][0], download['period_start'][1], download['period_start'][2]),\n",
    "                                     date(download['period_stop'][0], download['period_stop'][1], download['period_stop'][2])),\n",
    "                             platformname = 'Sentinel-1',\n",
    "                             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).sort_values('ingestiondate', ascending=[False])\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",
    "    # plot available Sentinel-1 tiles\n",
    "    geo_json = ipyleaflet.GeoJSON(data = products_json,\n",
    "                                  name = 'S1 tiles',\n",
    "                                  style = {'color' : 'royalblue'},\n",
    "                                  hover_style = {'fillOpacity' : 0.4})\n",
    "    download_map.add_layer(geo_json)\n",
    "    # call click_handler function when user clicks on Sentinel-1 tile\n",
    "    tile_ID = geo_json.on_hover(hover_handler)\n",
    "    # print table with download buttons\n",
    "    grid = ipywidgets.GridspecLayout(len(products_df.index)+1, 5)\n",
    "    grid[0,0] = ipywidgets.HTML('<h4>Index</h4>')\n",
    "    grid[0,1] = ipywidgets.HTML('<h4>Date</h4>')\n",
    "    grid[0,2] = ipywidgets.HTML('<h4>Polarisation</h4>')\n",
    "    grid[0,3] = ipywidgets.HTML('<h4>Size</h4>')\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",
    "# show product index and date on HTML element when hovering on Sentinel-1 tile \n",
    "def hover_handler(feature, **kwargs):\n",
    "    # reset HTML widget and set new values\n",
    "    html.value = '''\n",
    "    Index: {}<br/>\n",
    "    <small>Date: {}</small>\n",
    "    '''.format(feature['properties']['index'], feature['properties']['beginposition'])\n",
    "\n",
    "# search for available Sentinel-1 tiles with regard to AOI\n",
    "def on_searchButton_clicked(b):\n",
    "    # refocus of map to cover Sentinel-1 tiles\n",
    "    download_map.zoom = 6.5\n",
    "    # create WKT element with geographic coordinates of AOI\n",
    "    coordinates = draw_control.last_draw['geometry']['coordinates'][0]\n",
    "    lng1, lat1 = coordinates[0][0], coordinates[0][1]\n",
    "    lng2, lat2 = coordinates[1][0], coordinates[1][1]\n",
    "    lng3, lat3 = coordinates[2][0], coordinates[2][1]\n",
    "    lng4, lat4 = coordinates[3][0], coordinates[3][1]\n",
    "    footprint = 'POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))' % (lng1, lat1, lng2, lat2, lng3, lat3, lng4, lat4, lng1, lat1)\n",
    "    # call function to search for available Sentinel-1 tiles\n",
    "    queri(footprint)\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 download['imageDownload']:\n",
    "    # connect to the API\n",
    "    api = SentinelAPI(download['username'], download['password'], 'https://scihub.copernicus.eu/dhus')\n",
    "    # create map with search functionality\n",
    "    download_map = ipyleaflet.Map(zoom = 1.4)\n",
    "    download_map.add_control(ipyleaflet.SearchControl(\n",
    "        position = 'topleft',\n",
    "        url = 'https://nominatim.openstreetmap.org/search?format=json&q={s}',\n",
    "        zoom = 5))\n",
    "    download_map.add_control(ipyleaflet.ScaleControl(position='bottomleft'))\n",
    "    display(download_map)\n",
    "    \n",
    "    # create HTML widget to display product index and date of selected Sentinel-1 tile\n",
    "    html = ipywidgets.HTML('')\n",
    "    control = ipyleaflet.WidgetControl(widget=html, position='topright')\n",
    "    download_map.add_control(control)\n",
    "    \n",
    "    # check whether AOI file is given and convert to JSON in order to show AOI on map\n",
    "    try:\n",
    "        # calls readJSONfromAOI function to get GeoJSON from either JSON, SHP, or KMZ file\n",
    "        data_json = readJSONFromAOI('%s/AOI' % directory)\n",
    "        # show AOI on map according to JSON data\n",
    "        aoi = ipyleaflet.GeoJSON(data = data_json, style = {'color' : 'green'})\n",
    "        download_map.zoom = 6.5\n",
    "        try:\n",
    "            # GeoJSON format if KMZ is given\n",
    "            download_map.center = (aoi.data['features'][0]['geometry']['coordinates'][0][0][0][1],\n",
    "                                   aoi.data['features'][0]['geometry']['coordinates'][0][0][0][0])\n",
    "        except:\n",
    "            # GeoJSON format if JSON or SHP is given\n",
    "            download_map.center = (aoi.data['features'][0]['geometry']['coordinates'][0][0][1],\n",
    "                                   aoi.data['features'][0]['geometry']['coordinates'][0][0][0])\n",
    "        download_map.add_layer(aoi)\n",
    "        # search for available Sentinel-1 tiles according to JSON data\n",
    "        footprint = geojson_to_wkt(data_json)\n",
    "        queri(footprint)\n",
    "\n",
    "    # if no AOI is given, it needs to be selected manually\n",
    "    except FileNotFoundError:\n",
    "        # create draw control function\n",
    "        draw_control = ipyleaflet.DrawControl(circlemarker={},\n",
    "                                              polyline={},\n",
    "                                              polygon={'shapeOptions':{'color':'green'}},\n",
    "                                              rectangle={'shapeOptions':{'color':'green'}})\n",
    "        # adds draw control function to map\n",
    "        download_map.add_control(draw_control)\n",
    "        # create search button and call function to search for available Sentinel-1 tiles when clicked\n",
    "        searchButton = ipywidgets.Button(description = 'Search')\n",
    "        searchButton.on_click(on_searchButton_clicked)\n",
    "        display(searchButton)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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 drawing the area of interest. 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": {
    "code_folding": [
     0,
     7,
     61,
     98,
     118,
     129
    ]
   },
   "outputs": [],
   "source": [
    "# Click to run\n",
    "\n",
    "####################################################\n",
    "############### FUNCTION DEFINITIONS ###############\n",
    "####################################################\n",
    "\n",
    "# create S1 product\n",
    "def getScene(path):\n",
    "    # set correct path of input file and create S1 product\n",
    "    if len(files) is 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",
    "    center = (meta_data.getAttributeDouble('centre_lat'), meta_data.getAttributeDouble('centre_lon'))\n",
    "    locations = [[{'lat' : meta_data.getAttributeDouble('first_near_lat'), 'lng' : meta_data.getAttributeDouble('first_near_long')},\n",
    "                  {'lat' : meta_data.getAttributeDouble('last_near_lat'),  'lng' : meta_data.getAttributeDouble('last_near_long')},\n",
    "                  {'lat' : meta_data.getAttributeDouble('last_far_lat'),   'lng' : meta_data.getAttributeDouble('last_far_long')},\n",
    "                  {'lat' : meta_data.getAttributeDouble('first_far_lat'),  'lng' : meta_data.getAttributeDouble('first_far_long')}]]\n",
    "\n",
    "    # creates interactive map\n",
    "    basic_map = ipyleaflet.Map(center = center, zoom = 7.5)\n",
    "    # defines fixed polygon illustrating Sentinel-1 image\n",
    "    polygon_fix = ipyleaflet.Polygon(locations = locations, color='royalblue')\n",
    "    basic_map.add_layer(polygon_fix)\n",
    "    # displays map\n",
    "    basic_map.add_control(ipyleaflet.ScaleControl(position='bottomleft'))\n",
    "    display(basic_map)\n",
    "    \n",
    "    # check whether AOI file is given and convert to JSON in order to show AOI on map\n",
    "    try:\n",
    "        # calls readJSONfromAOI function to get GeoJSON from either JSON, SHP, or KMZ file\n",
    "        data_json = readJSONFromAOI('%s/AOI' % directory)\n",
    "        # show AOI on map according to JSON data\n",
    "        basic_map.add_layer(ipyleaflet.GeoJSON(data = data_json, style = {'color' : 'green'}))\n",
    "        # apply subset according to JSON data\n",
    "        footprint = geojson_to_wkt(data_json)\n",
    "        # run processing process\n",
    "        processing(S1_source, footprint)\n",
    "        \n",
    "    # if no AOI is given, it needs to be selected manually\n",
    "    except:\n",
    "        # editable polygon determining AOI\n",
    "        polygon_flex = ipyleaflet.Polygon(locations = locations,\n",
    "                                          color = 'green',\n",
    "                                          fill_color = 'green',\n",
    "                                          transform = True)\n",
    "        basic_map.add_layer(polygon_flex)\n",
    "        # create process button and call function to start ptocessing when clicked\n",
    "        processButton = ipywidgets.Button(description = 'Start Processing')\n",
    "        processButton.on_click(functools.partial(on_processButton_clicked,\n",
    "                                                 S1_source = S1_source,\n",
    "                                                 polygon_flex = polygon_flex))\n",
    "        display(processButton)\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",
    "\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 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",
    "    \n",
    "        numPixMinimum = len(band_data[abs(band_data - threshold_minimum) < 0.1])\n",
    "\n",
    "    # select 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",
    "# 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",
    "# start processing button \n",
    "def on_processButton_clicked(b, S1_source, polygon_flex):\n",
    "    # get coordinates and store in WKT format variable\n",
    "    lng1, lat1 = polygon_flex.locations[0][0]['lng'], polygon_flex.locations[0][0]['lat']\n",
    "    lng2, lat2 = polygon_flex.locations[0][1]['lng'], polygon_flex.locations[0][1]['lat']\n",
    "    lng3, lat3 = polygon_flex.locations[0][2]['lng'], polygon_flex.locations[0][2]['lat']\n",
    "    lng4, lat4 = polygon_flex.locations[0][3]['lng'], polygon_flex.locations[0][3]['lat']\n",
    "    footprint = 'POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))' % (lng1, lat1, lng2, lat2, lng3, lat3, lng4, lat4, lng1, lat1)\n",
    "    # run processing\n",
    "    processing(S1_source, footprint)\n",
    "\n",
    "# processing steps\n",
    "def processing(S1_source, footprint):\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",
    "    # 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 plotResoluts:\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 polarisations == 'both':\n",
    "    sourceBands = 'Amplitude_VH,Intensity_VH,Amplitude_VV,Intensity_VV'\n",
    "    output_extensions   = 'processed_VHVV'\n",
    "elif polarisations == 'VH':\n",
    "    sourceBands = 'Amplitude_VH,Intensity_VH'\n",
    "    output_extensions   = 'processed_VH'\n",
    "elif polarisations == '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",
    "    getScene('%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(getScene)\n",
    "    display(fc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Click to run\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) is not 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], polarisations)\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], polarisations))\n",
    "    ds = gdal.VectorTranslate('%s/%s_processed_%s.kml' % (KML_path, os.path.splitext(input_name)[0], polarisations), 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], polarisations))\n",
    "    shp_file.to_file('%s/%s_processed_%s.json' % (GeoJSON_path, os.path.splitext(input_name)[0], polarisations), 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",
    "results_map = ipyleaflet.Map(zoom=9, basemap=ipyleaflet.basemaps.OpenStreetMap.Mapnik)    \n",
    "if open_image.RasterCount == 1:\n",
    "    file = '%s/%s_processed_%s.json' % (GeoJSON_path, os.path.splitext(input_name)[0], polarisations)\n",
    "    with open(file, 'r') as f:\n",
    "        data_json = json.load(f) \n",
    "    mask = ipyleaflet.GeoJSON(data = data_json, name = 'Flood Mask', style = {'color':'blue', 'opacity':'1', 'fillColor':'blue', 'fillOpacity':'1', 'weight':'0.8'})\n",
    "    results_map.add_layer(mask)\n",
    "    results_map.center = (mask.data['features'][0]['geometry']['coordinates'][0][0][1],\n",
    "                          mask.data['features'][0]['geometry']['coordinates'][0][0][0])\n",
    "else:\n",
    "    file_VV = '%s/%s_processed_VV.json' % (GeoJSON_path, os.path.splitext(input_name)[0])\n",
    "    with open(file_VV, 'r') as f_VV:\n",
    "        data_json_VV = json.load(f_VV)\n",
    "    mask_VV = ipyleaflet.GeoJSON(data = data_json_VV, name = 'Flood Mask: VV', style = {'color':'red', 'opacity':'1', 'fillColor':'red', 'fillOpacity':'1', 'weight':'0.8'})\n",
    "    results_map.add_layer(mask_VV)\n",
    "    results_map.center = (mask_VV.data['features'][0]['geometry']['coordinates'][0][0][1],\n",
    "                          mask_VV.data['features'][0]['geometry']['coordinates'][0][0][0])  \n",
    "    file_VH = '%s/%s_processed_VH.json' % (GeoJSON_path, os.path.splitext(input_name)[0])\n",
    "    with open(file_VH, 'r') as f_VH:\n",
    "        data_json_VH = json.load(f_VH)\n",
    "    mask_VH = ipyleaflet.GeoJSON(data = data_json_VH, name = 'Flood Mask: VH', style = {'color':'blue', 'opacity':'1', 'fillColor':'blue', 'fillOpacity':'1', 'weight':'0.8'})\n",
    "    results_map.add_layer(mask_VH)\n",
    "results_map.add_control(ipyleaflet.FullScreenControl())\n",
    "results_map.add_control(ipyleaflet.LayersControl(position='topright'))\n",
    "results_map.add_control(ipyleaflet.ScaleControl(position='bottomleft'))\n",
    "display(results_map)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}