HYDRO30Workflow-v2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://radar.community.uaf.edu/wp-content/uploads/sites/667/2021/03/HydroSARbanner.jpg\" width=\"100%\" />\n",
    "<font face=\"Calibri\">\n",
    "<br>\n",
    "<font size=\"6\"> <b>HYDRO30 Product: Flood Mapping from Single Sentinel-1 SAR Images</b><img style=\"padding: 7px\" src=\"https://radar.community.uaf.edu/wp-content/uploads/sites/667/2021/03/UAFLogo_A_647.png\" width=\"170\" align=\"right\"/></font>\n",
    "\n",
    "<br>\n",
    "<font size=\"4\"> <b> Franz J Meyer, University of Alaska Fairbanks</b> <br>\n",
    "</font>\n",
    "\n",
    "<font size=\"3\">This approach is based on a methodology developed by the German Aerospace Center and published in <cite><a href=\"https://www.tandfonline.com/doi/full/10.1080/01431161.2016.1192304\"><i>Sentinel-1-based flood mapping: a fully automated processing chain</i></a> by Twele et al.</cite>. It is based on single-image Sentinel-1 SAR data and applies a dynamic thresholding method followed by fuzzy-logic-based post processing procedure. Modifications were made to the original approach to improve robustness and correct issues where the approach tended to fail.\n",
    "</font>\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# General Methodology and Workflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"3\">The workflow of the Sentinel-1-based processing chain, as outlined in the figure below, is composed of the following main elements: \n",
    "- <b>Step 1:</b> Automatic fetching of newly incoming SAR data over our area of interest and automatic processing to radiometric terrain corrected (RTC) products using the <a href=\"https://hyp3.asf.alaska.edu/\">ASF HyP3 service</a>.  \n",
    "- <b>Step 2:</b> initial classification using automatic thresholding, resulting in an initial flood mapping product (<b>Step 3</b>).\n",
    "- <b>Step 4:</b> fuzzy-logic-based classification refinement, \n",
    "- <b>Step 5:</b> final classification including auxiliary data, and \n",
    "- <b>Step 6:</b> dissemination of the results. \n",
    "\n",
    "\n",
    "<img style=\"padding: 7px\" src=\"https://courses.edx.org/asset-v1:AlaskaX+SAR-401+3T2020+type@asset+block@Watermappingworkflow2.jpg\" width=\"80%\" align=\"center\"/>\n",
    "</font>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flood Mapping Procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Loading Python Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Loading a number of Phython Libraries that are used in this Notebook\n",
    "import os\n",
    "import glob\n",
    "import osr\n",
    "import json\n",
    "import scipy\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "\n",
    "from osgeo import gdal\n",
    "from typing import Tuple\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plb # for add_patch, add_subplot, figure, hist, imshow, set_title, xaxis,_label, text \n",
    "import matplotlib.pyplot as plt # for add_subplot, axis, figure, imshow, legend, plot, set_axis_off, set_data,\n",
    "                                # set_title, set_xlabel, set_ylabel, set_ylim, subplots, title, twinx\n",
    "import matplotlib.patches as patches  # for Rectangle\n",
    "import matplotlib.animation as an # for FuncAnimation\n",
    "from matplotlib import rc \n",
    "\n",
    "import ipywidgets as ui\n",
    "import pylab as pl\n",
    "\n",
    "try:\n",
    "    import rasterio\n",
    "except:\n",
    "    !pip install rasterio\n",
    "    import rasterio\n",
    "try:\n",
    "    import pyproj\n",
    "except:\n",
    "    !pip install pyproj\n",
    "    import pyproj\n",
    "    \n",
    "try:\n",
    "    import skfuzzy\n",
    "except:   \n",
    "    !pip install scikit-fuzzy\n",
    "    import skfuzzy\n",
    "    \n",
    "try:\n",
    "    from tqdm.auto import tqdm\n",
    "except:\n",
    "    !pip install tqdm\n",
    "    from tqdm.auto import tqdm\n",
    "\n",
    "try:\n",
    "    import astropy\n",
    "except:\n",
    "    !pip install --user astropy\n",
    "    import astropy\n",
    "import astropy.convolution\n",
    "from scipy import ndimage\n",
    "\n",
    "import warnings #Suppress warnings on occasion\n",
    "    \n",
    "from asf_notebook import new_directory\n",
    "from asf_notebook import handle_old_data\n",
    "from asf_notebook import path_exists\n",
    "from asf_notebook import input_path\n",
    "\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Helper Scripts To Set Up the Flood Mapping Approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function to pad an image, so it may be split into tiles with consistent dimensions</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Function to Pad the image\n",
    "def pad_image(image: np.ndarray, to: int) -> np.ndarray:\n",
    "    height, width = image.shape\n",
    "\n",
    "    n_rows, n_cols = get_tile_row_col_count(height, width, to)\n",
    "    new_height = n_rows * to\n",
    "    new_width = n_cols * to\n",
    "\n",
    "    padded = np.zeros((new_height, new_width))\n",
    "    padded[:image.shape[0], :image.shape[1]] = image\n",
    "    return padded"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function to tile an image</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Function for Image Tiling\n",
    "def tile_image(image: np.ndarray, width, height) -> np.ndarray:\n",
    "    _nrows, _ncols = image.shape\n",
    "    _strides = image.strides\n",
    "\n",
    "    nrows, _m = divmod(_nrows, height)\n",
    "    ncols, _n = divmod(_ncols, width)\n",
    "\n",
    "    assert _m == 0, \"Image must be evenly tileable. Please pad it first\"\n",
    "    assert _n == 0, \"Image must be evenly tileable. Please pad it first\"\n",
    "\n",
    "    return np.lib.stride_tricks.as_strided(\n",
    "        np.ravel(image),\n",
    "        shape=(nrows, ncols, height, width),\n",
    "        strides=(height * _strides[0], width * _strides[1], *_strides),\n",
    "        writeable=False\n",
    "    ).reshape(nrows * ncols, height, width)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function for multi-class Expectation Maximization Thresholding.</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Function for Threshold Calculation using an Expectation Maximization Approach\n",
    "def EMSeg_opt(image, number_of_classes):\n",
    "    image_copy = image.copy()\n",
    "    image_copy2 = np.ma.filled(image.astype(float), np.nan) # needed for valid posterior_lookup keys\n",
    "    image = image.flatten()\n",
    "    minimum = np.amin(image)\n",
    "    image = image - minimum + 1\n",
    "    maximum = np.amax(image)\n",
    "\n",
    "    size = image.size\n",
    "    histogram = make_histogram(image)\n",
    "    nonzero_indices = np.nonzero(histogram)[0]\n",
    "    histogram = histogram[nonzero_indices]\n",
    "    histogram = histogram.flatten()\n",
    "    class_means = (\n",
    "            (np.arange(number_of_classes) + 1) * maximum /\n",
    "            (number_of_classes + 1)\n",
    "    )\n",
    "    class_variances = np.ones((number_of_classes)) * maximum\n",
    "    class_proportions = np.ones((number_of_classes)) * 1 / number_of_classes\n",
    "    sml = np.mean(np.diff(nonzero_indices)) / 1000\n",
    "    iteration = 0\n",
    "    while(True):\n",
    "        class_likelihood = make_distribution(\n",
    "            class_means, class_variances, class_proportions, nonzero_indices\n",
    "        )\n",
    "        sum_likelihood = np.sum(class_likelihood, 1) + np.finfo(\n",
    "                class_likelihood[0][0]).eps\n",
    "        log_likelihood = np.sum(histogram * np.log(sum_likelihood))\n",
    "        for j in range(0, number_of_classes):\n",
    "            class_posterior_probability = (\n",
    "                histogram * class_likelihood[:,j] / sum_likelihood\n",
    "            )\n",
    "            class_proportions[j] = np.sum(class_posterior_probability)\n",
    "            class_means[j] = (\n",
    "                np.sum(nonzero_indices * class_posterior_probability)\n",
    "                    / class_proportions[j]\n",
    "            )\n",
    "            vr = (nonzero_indices - class_means[j])\n",
    "            class_variances[j] = (\n",
    "                np.sum(vr *vr * class_posterior_probability)\n",
    "                    / class_proportions[j] +sml\n",
    "            )\n",
    "            del class_posterior_probability, vr\n",
    "        class_proportions = class_proportions + 1e-3\n",
    "        class_proportions = class_proportions / np.sum(class_proportions)\n",
    "        class_likelihood = make_distribution(\n",
    "            class_means, class_variances, class_proportions, nonzero_indices\n",
    "        )\n",
    "        sum_likelihood = np.sum(class_likelihood, 1) + np.finfo(\n",
    "                class_likelihood[0,0]).eps\n",
    "        del class_likelihood\n",
    "        new_log_likelihood = np.sum(histogram * np.log(sum_likelihood))\n",
    "        del sum_likelihood\n",
    "        if((new_log_likelihood - log_likelihood) < 0.000001):\n",
    "            break\n",
    "        iteration = iteration + 1\n",
    "    del log_likelihood, new_log_likelihood\n",
    "    class_means = class_means + minimum - 1\n",
    "    s = image_copy.shape\n",
    "    posterior = np.zeros((s[0], s[1], number_of_classes))\n",
    "    posterior_lookup = dict()\n",
    "    for i in range(0, s[0]):\n",
    "        for j in range(0, s[1]):\n",
    "            pixel_val = image_copy2[i,j] \n",
    "            if pixel_val in posterior_lookup:\n",
    "                for n in range(0, number_of_classes): \n",
    "                    posterior[i,j,n] = posterior_lookup[pixel_val][n]\n",
    "            else:\n",
    "                posterior_lookup.update({pixel_val: [0]*number_of_classes})\n",
    "                for n in range(0, number_of_classes): \n",
    "                    x = make_distribution(\n",
    "                        class_means[n], class_variances[n], class_proportions[n],\n",
    "                        image_copy[i,j]\n",
    "                    )\n",
    "                    posterior[i,j,n] = x * class_proportions[n]\n",
    "                    posterior_lookup[pixel_val][n] = posterior[i,j,n]\n",
    "    return posterior, class_means, class_variances, class_proportions\n",
    "\n",
    "def make_histogram(image):\n",
    "    image = image.flatten()\n",
    "    indices = np.nonzero(np.isnan(image))\n",
    "    image[indices] = 0\n",
    "    indices = np.nonzero(np.isinf(image))\n",
    "    image[indices] = 0\n",
    "    del indices\n",
    "    size = image.size\n",
    "    maximum = int(np.ceil(np.amax(image)) + 1)\n",
    "    #maximum = (np.ceil(np.amax(image)) + 1)\n",
    "    histogram = np.zeros((1, maximum))\n",
    "    for i in range(0,size):\n",
    "        #floor_value = int(np.floor(image[i]))\n",
    "        floor_value = np.floor(image[i]).astype(np.uint8)\n",
    "        #floor_value = (np.floor(image[i]))\n",
    "        if floor_value > 0 and floor_value < maximum - 1:\n",
    "            temp1 = image[i] - floor_value\n",
    "            temp2 = 1 - temp1\n",
    "            histogram[0,floor_value] = histogram[0,floor_value] + temp1\n",
    "            histogram[0,floor_value - 1] = histogram[0,floor_value - 1] + temp2\n",
    "    histogram = np.convolve(histogram[0], [1,2,3,2,1])\n",
    "    histogram = histogram[2:(histogram.size - 3)]\n",
    "    histogram = histogram / np.sum(histogram)\n",
    "    return histogram\n",
    "\n",
    "def make_distribution(m, v, g, x):\n",
    "    x = x.flatten()\n",
    "    m = m.flatten()\n",
    "    v = v.flatten()\n",
    "    g = g.flatten()\n",
    "    y = np.zeros((len(x), m.shape[0]))\n",
    "    for i in range(0,m.shape[0]):\n",
    "        d = x - m[i]\n",
    "        amp = g[i] / np.sqrt(2*np.pi*v[i])\n",
    "        y[:,i] = amp * np.exp(-0.5 * (d * d) / v[i])\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function to calculate the number of rows and columns of tiles needed to tile an image to a given size</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Function to Calculate Number of Image Tiles\n",
    "def get_tile_row_col_count(height: int, width: int, tile_size: int) -> Tuple[int, int]:\n",
    "    return int(np.ceil(height / tile_size)), int(np.ceil(width / tile_size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function to extract the tiff dates from a wildcard path:</b> </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Function to Extract Image Acquisition Date Information\n",
    "def get_dates(paths):\n",
    "    dates = []\n",
    "    pths = glob.glob(paths)\n",
    "    for p in pths:\n",
    "        date = p.split('/')[-1].split(\"_\")[3].split(\"T\")[0]\n",
    "        dates.append(date)\n",
    "    dates.sort()\n",
    "    return dates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function to save a mask</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Some Additional Python Utilities\n",
    "def write_mask_to_file(mask: np.ndarray, file_name: str, projection: str, geo_transform: str) -> None:\n",
    "    (width, height) = mask.shape\n",
    "    out_image = gdal.GetDriverByName('GTiff').Create(\n",
    "        file_name, height, width, bands=1\n",
    "    )\n",
    "    out_image.SetProjection(projection)\n",
    "    out_image.SetGeoTransform(geo_transform)\n",
    "    out_image.GetRasterBand(1).WriteArray(mask)\n",
    "    out_image.GetRasterBand(1).SetNoDataValue(0)\n",
    "    out_image.FlushCache()\n",
    "\n",
    "def gdal_write(ary, geoTransform, fileformat=\"GTiff\", filename='jupyter_rocks.tif', format=gdal.GDT_Float64, nodata=None, srs_proj4='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'):\n",
    "    '''gdal_write(ary, geoTransform, format=\"GTiff\", filename='jupyter_rocks.tif', format=gdal.GDT_Float64 nodata=None, srs_proj4='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')\n",
    "    ary: 2D array.\n",
    "    geoTransform: [top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution]\n",
    "    format: \"GTiff\"     \n",
    "    '''           \n",
    "    if ary.ndim ==2:\n",
    "      Ny, Nx = ary.shape\n",
    "      Nb = 1;\n",
    "    elif ary.ndim==3:\n",
    "      Ny,Nx,Nb=ary.shape\n",
    "    else: \n",
    "      print(\"Input array has to be 2D or 3D.\")\n",
    "      return None\n",
    "    \n",
    "    driver = gdal.GetDriverByName(fileformat)\n",
    "    ds = driver.Create(filename, Nx, Ny, Nb, gdal.GDT_Float64)\n",
    "\n",
    "    #ds.SetGeoTransform( ... ) # define GeoTransform tuple\n",
    "    # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution\n",
    "    ds.SetGeoTransform( geoTransform )    \n",
    "    srs=osr.SpatialReference()\n",
    "    srs.ImportFromProj4(srs_proj4)\n",
    "    ds.SetProjection(srs.ExportToWkt() );\n",
    "    if nodata is not None:\n",
    "        ds.GetRasterBand(1).SetNoDataValue(0);\n",
    "    if Nb==1:\n",
    "      ds.GetRasterBand(1).WriteArray(ary)\n",
    "    else:\n",
    "      for b in range(Nb):\n",
    "        ds.GetRasterBand(b+1).WriteArray(ary[:,:,b])\n",
    "    ds = None\n",
    "    #print(\"File written to: \" + filename);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Function To Select Path to Folder Holding your HAND File\n",
    "class PathSelector():\n",
    "\n",
    "    def __init__(self,start_dir,select_file=True):\n",
    "        self.file        = None \n",
    "        self.select_file = select_file\n",
    "        self.cwd         = start_dir\n",
    "        self.select      = ui.SelectMultiple(options=['init'],value=(),rows=10,description='') \n",
    "        self.accord      = ui.Accordion(children=[self.select]) \n",
    "\n",
    "        self.accord.selected_index = None # Start closed (showing path only)\n",
    "        self.refresh(self.cwd)\n",
    "        self.select.observe(self.on_update,'value')\n",
    "\n",
    "    def on_update(self,change):\n",
    "        if len(change['new']) > 0:\n",
    "            self.refresh(change['new'][0])\n",
    "\n",
    "    def refresh(self,item):\n",
    "        path = os.path.abspath(os.path.join(self.cwd,item))\n",
    "\n",
    "        if os.path.isfile(path):\n",
    "            if self.select_file:\n",
    "                self.accord.set_title(0,path)  \n",
    "                self.file = path\n",
    "                self.accord.selected_index = None\n",
    "            else:\n",
    "                self.select.value = ()\n",
    "\n",
    "        else: # os.path.isdir(path)\n",
    "            self.file = None \n",
    "            self.cwd  = path\n",
    "\n",
    "            # Build list of files and dirs\n",
    "            keys = ['[..]']; \n",
    "            for item in os.listdir(path):\n",
    "                if item[0] == '.':\n",
    "                    continue\n",
    "                elif os.path.isdir(os.path.join(path,item)):\n",
    "                    keys.append('['+item+']'); \n",
    "                else:\n",
    "                    keys.append(item); \n",
    "\n",
    "            # Sort and create list of output values\n",
    "            keys.sort(key=str.lower)\n",
    "            vals = []\n",
    "            for k in keys:\n",
    "                if k[0] == '[':\n",
    "                    vals.append(k[1:-1]) # strip off brackets\n",
    "                else:\n",
    "                    vals.append(k)\n",
    "\n",
    "            # Update widget\n",
    "            self.accord.set_title(0,path)  \n",
    "            self.select.options = list(zip(keys,vals)) \n",
    "            with self.select.hold_trait_notifications():\n",
    "                self.select.value = ()\n",
    "                \n",
    "def get_proj4(filename):\n",
    "    f=rasterio.open(filename)\n",
    "    return pyproj.Proj(f.crs, preserve_units=True)  #used in pysheds\n",
    "    \n",
    "def gdal_read(filename, ndtype=np.float64):\n",
    "    '''\n",
    "    z=readData('/path/to/file')\n",
    "    '''\n",
    "    ds = gdal.Open(filename) \n",
    "    return np.array(ds.GetRasterBand(1).ReadAsArray()).astype(ndtype);\n",
    "def gdal_get_geotransform(filename):\n",
    "    '''\n",
    "    [top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution]=gdal_get_geotransform('/path/to/file')\n",
    "    '''\n",
    "    #http://stackoverflow.com/questions/2922532/obtain-latitude-and-longitude-from-a-geotiff-file\n",
    "    ds = gdal.Open(filename)\n",
    "    return ds.GetGeoTransform()\n",
    "\n",
    "def fill_nan(arr):\n",
    "    \"\"\"\n",
    "    filled_arr=fill_nan(arr)\n",
    "    Fills Not-a-number values in arr using astropy. \n",
    "    \"\"\"    \n",
    "    kernel = astropy.convolution.Gaussian2DKernel(x_stddev=3) #kernel x_size=8*stddev\n",
    "    arr_type=arr.dtype          \n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter(\"ignore\")\n",
    "        while np.any(np.isnan(arr)):\n",
    "            arr = astropy.convolution.interpolate_replace_nans(arr.astype(float), kernel, convolve=astropy.convolution.convolve)\n",
    "    return arr.astype(arr_type) \n",
    "\n",
    "def get_tiff_paths(paths: str) -> list:\n",
    "    tiff_paths = !ls $paths | sort -t_ -k5,5\n",
    "    return tiff_paths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Find SAR Data Sets to Process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Enter the path to the data stack</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Define Path to Folder Containing Your SAR Images\n",
    "while True:\n",
    "    print(\"Enter the absolute path to the directory holding your tiffs.\")\n",
    "    tiff_dir = input()\n",
    "    paths = f\"{tiff_dir}/*_V*.tif*\"\n",
    "    if os.path.exists(tiff_dir):\n",
    "        tiff_paths = get_tiff_paths(paths)\n",
    "        if len(tiff_paths) < 1:\n",
    "            print(f\"{tiff_dir} exists but contains no tifs.\")\n",
    "            print(\"You will not be able to proceed until tifs are prepared.\")\n",
    "        break\n",
    "    else:\n",
    "        print(f\"\\n{tiff_dir} does not exist.\")\n",
    "        continue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Move into the parent directory of the directory containing the data and create a directory in which to store the water masks</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Move into Data Directory for Processing\n",
    "analysis_directory = os.path.dirname(tiff_dir)\n",
    "os.chdir(analysis_directory)\n",
    "mask_directory = f'{analysis_directory}/Water_Masks'\n",
    "new_directory(mask_directory)\n",
    "print(f\"Current working directory: {os.getcwd()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function to create a dictionary containing lists of each vv/vh pair</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Define Function for Grouping VV / VH \n",
    "def group_polarizations(tiff_paths: list) -> dict:\n",
    "    pths = {}\n",
    "    for tiff in tiff_paths:\n",
    "        product_name = tiff.split('.')[0][:-2]\n",
    "        if product_name in pths:\n",
    "            pths[product_name].append(tiff)\n",
    "        else:\n",
    "            pths.update({product_name: [tiff]})\n",
    "            pths[product_name].sort()\n",
    "    return pths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Write a function to confirm the presence of both VV and VH images in all image sets</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Define Function to Confirm that every Image Date has Both VV and VH Polarization\n",
    "def confirm_dual_polarizations(paths: dict) -> bool:\n",
    "    for p in paths:\n",
    "        if len(paths[p]) == 2:\n",
    "            if ('vv' not in paths[p][1] and 'VV' not in paths[p][1]) or \\\n",
    "            ('vh' not in paths[p][0] and 'VH' not in paths[p][0]):\n",
    "                return False\n",
    "    return True   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><b>Create a dictionary of VV/VH pairs and check it for completeness</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Group VV and VH Images by Date and make sure Every Date has both Polarizations \n",
    "grouped_pths = group_polarizations(tiff_paths)\n",
    "if not confirm_dual_polarizations(grouped_pths):\n",
    "    print(\"ERROR: AI_Water requires both VV and VH polarizations.\")\n",
    "else:\n",
    "    print(\"Confirmed presence of VV and VH polarities for each product.\")\n",
    "    \n",
    "#print(grouped_pths) #uncomment to print VV/VH path pairs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Identify HAND Layer to be used in Flood Mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"><font color='rgba(200,0,0,0.2)'><b>Note:</b></font> Please pick a HAND layer that covers your entire area of interest (or more).</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Pick the HAND File you want to use\n",
    "print(\"Choose your HAND layer:\")\n",
    "f = PathSelector('.')\n",
    "display(f.accord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Read HAND Data and Projection Information\n",
    "HAND_file=f.accord.get_title(0)\n",
    "info = (gdal.Info(HAND_file, options = ['-json']))\n",
    "info = json.dumps(info)\n",
    "info = (json.loads(info))['coordinateSystem']['wkt']\n",
    "Hzone = info.split('ID')[-1].split(',')[1][0:-2]\n",
    "Hproj_type = info.split('ID')[-1].split('\"')[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Do Initial Flood Mapping using Adaptive Dynamic Thresholding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"4\"><b>A bit of background on the approach:</b></font>\n",
    "\n",
    "<font face=\"Calibri\" size=\"3\"><u><b>This is what's implemented so far to arrive at the initial flood mapping product </b> (approach modified from Twele et al., (2016))</u>: An automatic tile-based thresholding procedure (Martinis, Twele, and Voigt 2009) is used to generate an initial land/water classification. The selection of tiles is performed on a bilevel quadtree structure with parent level $L^+$ and child level $L^−$:\n",
    "    \n",
    "1. Firstly, the entire data are separated into quadratic non-overlapping parent tiles on level $L^+$ with a size of $200 \\times 200$ pixels. Each parent object is further represented by four quadratic child objects on a second level $L^−$. The tile selection process is based on statistical hierarchical relations between parent and child objects. <font color='rgba(200,0,0,0.2)'><b>[Note: I modified this step slightly. Instead of starting with a $200 \\times 200$ pixels tile, I am using $100 \\times 100$ pixel tiles instead. This seems to be be more robust. This could be modified in the future</b></font>\n",
    "2. A number of parent tiles is automatically selected which offer the highest (>95% quantile) standard deviation on $L^+$ of the mean backscatter values of the respective child objects on $L^−$. This criterion serves as a measure of the degree of variation within the data and can therefore be used as an indicator of the probability that the tiles are characterized by spatial inhomogeneity and contain more than one semantic class. <font color='rgba(200,0,0,0.2)'><b>[Note: I modified this step slightly. Instead of the standard deviation $\\sigma$, I am using the coefficient of variation $\\frac{\\sigma}{\\mu}$ here. This seems to be more stable.]</b></font> The selected parent objects should also have a mean individual backscatter value lower than the mean of all parent tiles on $L^+$. This ensures that tiles lying on the boundary between water and no water areas are selected. In case that no tiles fulfil these criteria, the tile size on $L^+$ and $L^−$ is halved and the quantile for the tile selection is reduced to 90% to guarantee a successful tile selection also in data with a relatively low extent of water surfaces or with smaller dispersed water bodies. \n",
    "3. <b>NEW - Integration of HAND</b>: To improve the robustness of the automatic threshold derivation it should be considered to restrict the tile selection in Step (3) to only pixels situated in flood-prone regions defined by a HAND-based binary Exclustion Mask (HAND-EM). To create HAND-EM, a threshold is applied to HAND to identify non-flood prone areas. A threshold value of $\\geq 15m$ is proposed. The HAND-EM further is shrunk by one pixel using an 8-neighbour function to account for potential geometric inaccuracies between the exclude layer and SAR data. Tiles are only considered in case less than 20% of its data pixels are excluded by HAND-EM. <font color='rgba(200,0,0,0.2)'><b>This step is now included!</b></font>\n",
    "4. Out of the number of the initially selected tiles, a limited number of N parent tiles are finally chosen for threshold computation. This selection is accomplished by ranking the parent tiles according to the standard deviation of the mean backscatter values of the respective child objects. Tiles with the highest values are chosen for $N$. Extensive testing yielded that $N = 5$ is a sufficient number of parent tiles for threshold computation. \n",
    "5. The parametric Kittler and Illingworth minimum error thresholding approach (Kittler and Illingworth 1986) is then employed to derive local threshold values using a cost function which is based on the statistical parameterization of the sub-histograms of all selected tiles as bi-modal Gaussian mixture distributions. In order to derive a global (i.e. scenebased) threshold, the locally derived thresholds are combined by computing their arithmetic mean. <font color='rgba(200,0,0,0.2)'><b>[Note: I replaced Kittler Illingworth with a multi-mode Expectation Maximization approach to be able to cope with regions that have three classes (urban and bright slopes; average brightness farm lands; water). This happens a lot in the midwest and Kittler Illingworth didn't perform well.]</b></font>\n",
    "6. Using the dynamically calculated threshold, both the VV and VH scenes are thresholded for water detection\n",
    "7. The detected water maps are combined to arrive at an intial water mask (<b>Step 3</b>) that is further refined in post processing (Section 2.6).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"4\"><b>Now Let's do the Work:</b></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     18
    ],
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Perform Inital Surface Water Mask Creation\n",
    "\n",
    "# define some variables you might want to change\n",
    "precentile = 0.95        # Standard deviation percentile threshold for pivotal tile selection\n",
    "tilesize = 100           # Setting tile size to use in thresholding algorithm\n",
    "tilesize2 = 50\n",
    "Hpick = 0.8              # Threshold for required fraction of valid HAND-EM pixels per tile\n",
    "vv_corr = -15.5          # VV threshold to use if threshold calculation did not succeed\n",
    "vh_corr = -23.0          # VH threshold to use if threshold calculation did not succeed\n",
    "Hthresh = 15\n",
    "\n",
    "# Now do adaptive threshold selection\n",
    "vv_thresholds = np.array([])\n",
    "vh_thresholds = np.array([])\n",
    "floodarea = np.array([])\n",
    "vh_thresholds_corr = np.array([])\n",
    "vv_thresholds_corr = np.array([])\n",
    "\n",
    "for pair in grouped_pths:\n",
    "    for tiff in grouped_pths[pair]:\n",
    "        resampled_dem_path=f'{tiff_dir}/resamp_dem.tif'\n",
    "        f = gdal.Open(tiff)\n",
    "        img_array = f.ReadAsArray()\n",
    "        original_shape = img_array.shape\n",
    "        img_array = 0     # free up RAM\n",
    "        \n",
    "        info1 = (gdal.Info(tiff, options = ['-json']))\n",
    "        info1 = json.dumps(info1)\n",
    "        ul = (json.loads(info1))['cornerCoordinates']['upperLeft']\n",
    "        lr = (json.loads(info1))['cornerCoordinates']['lowerRight']\n",
    "        coordsys = (json.loads(info1))['coordinateSystem']['wkt']\n",
    "        Szone = coordsys.split('ID')[-1].split(',')[1][0:-2]\n",
    "        Sproj_type = coordsys.split('ID')[-1].split('\"')[1]\n",
    "\n",
    "        west= ul[0]\n",
    "        east= lr[0]\n",
    "        south= lr[1]\n",
    "        north= ul[1]\n",
    "        \n",
    "        print('----------------------------------')\n",
    "        print('Extracting Relevant Subset from HAND Layer')\n",
    "        print('----------------------------------')\n",
    "        \n",
    "        cmd_resamp=f\"gdalwarp -overwrite -s_srs {Hproj_type}:\"\\\n",
    "        f\"{Hzone} -t_srs EPSG:{Szone} -te {west} {south} {east} {north} -ts {original_shape[1]} {original_shape[0]} -r lanczos {HAND_file} {resampled_dem_path}\"\n",
    "        #print(cmd_resamp)\n",
    "        os.system(cmd_resamp)\n",
    "        \n",
    "        g = gdal.Open(resampled_dem_path)\n",
    "        hand = g.ReadAsArray()\n",
    "        Hmask = hand < Hthresh\n",
    "        handem = np.zeros_like(hand)\n",
    "        sel = np.ones_like(hand)\n",
    "        handem[Hmask] = sel[Hmask]\n",
    "        hand = 0\n",
    "        \n",
    "        # Tile up HAND-EM data\n",
    "        handem_p = pad_image(handem, tilesize)\n",
    "        hand_tiles = tile_image(handem_p,width=tilesize,height=tilesize)\n",
    "        Hsum = np.sum(hand_tiles, axis=(1,2))\n",
    "        Hpercent = Hsum/(tilesize*tilesize)\n",
    "        \n",
    "        n_rows, n_cols = get_tile_row_col_count(*original_shape, tile_size=tilesize)\n",
    "        print('----------------------------------')\n",
    "        print(f'Flood Mapping on Image: {tiff}')\n",
    "        print('----------------------------------')\n",
    "        if 'vv' in tiff or 'VV' in tiff:\n",
    "            vv_array = pad_image(f.ReadAsArray(), tilesize)\n",
    "            invalid_pixels = np.nonzero(vv_array == 0.0)\n",
    "            vv_tiles = tile_image(vv_array,width=tilesize,height=tilesize)\n",
    "            a = np.shape(vv_tiles)\n",
    "            vv_std = np.zeros(a[0])\n",
    "            vvt_masked = np.ma.masked_where(vv_tiles<=0, vv_tiles)\n",
    "            vv_picktiles = np.zeros_like(vv_tiles)\n",
    "            for k in range(a[0]):\n",
    "                vv_subtiles = tile_image(vvt_masked[k,:,:],width=tilesize2,height=tilesize2)\n",
    "                vvst_mean = np.ma.mean(vv_subtiles, axis=(1,2))\n",
    "                vvst_std = np.ma.std(vvst_mean)\n",
    "                vv_std[k] = np.ma.std(vvst_mean) \n",
    "            \n",
    "            # find tiles with largest standard deviations\n",
    "            vv_mean = np.ma.median(vvt_masked, axis=(1,2))\n",
    "            x_vv = np.sort(vv_std/vv_mean)\n",
    "            y_vv = np.arange(1, x_vv.size+1) / x_vv.size\n",
    "            \n",
    "            percentile2 = precentile\n",
    "            noconverge_vv = 0\n",
    "            sort_index = 0\n",
    "            while np.size(sort_index) < 5 and noconverge_vv == 0: \n",
    "                threshold_index_vv = np.ma.min(np.where(y_vv>percentile2))\n",
    "                threshold_vv = x_vv[threshold_index_vv]\n",
    "                #sd_select_vv = np.nonzero(vv_std/vv_mean>threshold_vv)\n",
    "                s_select_vv = np.nonzero(vv_std/vv_mean>threshold_vv) \n",
    "                h_select_vv = np.nonzero(Hpercent > Hpick)               # Includes HAND-EM in selection\n",
    "                sd_select_vv = np.intersect1d(s_select_vv, h_select_vv)\n",
    "            \n",
    "                # find tiles with mean values lower than the average mean\n",
    "                omean_vv = np.ma.median(vv_mean[h_select_vv])\n",
    "                mean_select_vv = np.nonzero(vv_mean<omean_vv)\n",
    "            \n",
    "                # Intersect tiles with large std with tiles that have small means\n",
    "                msdselect_vv = np.intersect1d(sd_select_vv, mean_select_vv)\n",
    "                sort_index = np.flipud(np.argsort(vv_std[msdselect_vv]))\n",
    "                percentile2 = percentile2 - 0.01\n",
    "                if percentile2 < 0.5:\n",
    "                    noconverge_vv = 1\n",
    "                    sort_index = threshold_index_vv\n",
    "                    m_thresh_vv = (vv_corr/10.0+30)\n",
    "                    print(\"Tile Selection Did Not Converge - Going with Default Threshold\")\n",
    "            \n",
    "            if noconverge_vv < 1:\n",
    "                finalselect_vv = sort_index[0:5]\n",
    "            \n",
    "            temp = np.ma.masked_where(vv_array<=0, vv_array)\n",
    "            np.seterr(divide='ignore')\n",
    "            np.seterr(invalid='ignore')\n",
    "            dbvv = np.log10(temp)+30\n",
    "                \n",
    "            if noconverge_vv < 1:\n",
    "                # find local thresholds for 5 \"best\" tiles in the image\n",
    "                l_thresh_vv = np.zeros(5)\n",
    "                EMthresh_vv = np.zeros(5)\n",
    "                scaling = 256/(np.mean(dbvv) + 3*np.std(dbvv))\n",
    "                #scaling = 256/(np.mean(vv_array) + 3*np.std(vv_array))\n",
    "                np.seterr(divide='ignore')\n",
    "                np.seterr(invalid='ignore')\n",
    "                dbtile = np.log10(vvt_masked)+30\n",
    "                for k in range(5):\n",
    "                    test = dbtile[msdselect_vh[finalselect_vh[k]]] * scaling\n",
    "                    #test = vvt_masked[msdselect_vv[finalselect_vv[k]]] * scaling\n",
    "                    A = np.around(test)\n",
    "                    A = A.astype(int)\n",
    "                    #t_thresh = Kittler(A)\n",
    "                    [posterior, cm, cv, cp] = EMSeg_opt(A, 3)\n",
    "                    sorti = np.argsort(cm)\n",
    "                    cms = cm[sorti]\n",
    "                    cvs = cv[sorti]\n",
    "                    cps = cp[sorti]\n",
    "                    xvec = np.arange(cms[0],cms[1],step=.05)\n",
    "                    x1 = make_distribution(cms[0], cvs[0], cps[0], xvec)\n",
    "                    x2 = make_distribution(cms[1], cvs[1], cps[1], xvec)\n",
    "                    dx = np.abs(x1 - x2)\n",
    "                    diff1 = posterior[:,:,0] - posterior[:,:,1]\n",
    "                    t_ind = np.argmin(dx)\n",
    "                    EMthresh_vv[k] = xvec[t_ind]/scaling\n",
    "                \n",
    "                    #l_thresh_vv[k] = t_thresh / scaling\n",
    "                    #dbtile = np.log10(vvt_masked)+30\n",
    "                \n",
    "                    # Mark Tiles used for Threshold Estimation\n",
    "                    vv_picktiles[msdselect_vh[finalselect_vh[k]],:,:]= np.ones_like(vv_tiles[msdselect_vh[finalselect_vh[k]],:,:])\n",
    "            \n",
    "            if noconverge_vv < 1:\n",
    "                # Calculate best threshold for VV and VH as the mean of the 5 thresholds calculated in the previous section \n",
    "                #m_thresh_vv = np.median(l_thresh_vv)\n",
    "                #print(EMthresh_vv-30)\n",
    "                EMts = np.sort(EMthresh_vv)\n",
    "                #m_thresh_vv = np.median(EMthresh_vv)\n",
    "                m_thresh_vv = np.median(EMts[0:4])\n",
    "            \n",
    "            print(\"Best VV Flood Mapping Threshold [dB]: %.2f\" % (10*(m_thresh_vv-30)))\n",
    "            print(\" \")\n",
    "            \n",
    "            # Derive flood mask using the best threshold\n",
    "            if m_thresh_vv < (vv_corr/10.0+30):\n",
    "                change_mag_mask_vv = np.ma.masked_where(dbvv<=0, dbvv) < m_thresh_vv\n",
    "                vv_thresholds_corr = np.append(vv_thresholds_corr, 10.0*(m_thresh_vv-30))\n",
    "                #change_mag_mask_vv = np.ma.masked_where(vv_array==0, vv_array) < m_thresh_vv\n",
    "            else:\n",
    "                change_mag_mask_vv = np.ma.masked_where(dbvv<=0, dbvv) < (vv_corr/10.0+30)\n",
    "                vv_thresholds_corr = np.append(vv_thresholds_corr, vv_corr)\n",
    "            \n",
    "            # Create Binary masks showing flooded pixels as \"1\"s\n",
    "            flood_vv = np.zeros_like(vv_array)\n",
    "            sel = np.ones_like(vv_array)\n",
    "            flood_vv[change_mag_mask_vv] = sel[change_mag_mask_vv]\n",
    "            np.putmask(flood_vv,vv_array==0 , 0)\n",
    "            \n",
    "            # Export flood maps as GeoTIFFs\n",
    "            filename, ext = os.path.basename(tiff).split('.')\n",
    "            outfile = f\"{mask_directory}/{filename}_water_mask.{ext}\"\n",
    "            write_mask_to_file(flood_vv, outfile, f.GetProjection(), f.GetGeoTransform())\n",
    "            #vv_array = 0\n",
    "            maskedarray = 0\n",
    "            temp = 0\n",
    "            dbvv = 0\n",
    "            change_mag_mask_vv = 0\n",
    "            Hmask = 0\n",
    "            Hpercent = 0\n",
    "            Hsum = 0\n",
    "            hand_tiles = 0\n",
    "            handem = 0\n",
    "            handem_p = 0\n",
    "            vv_tiles = 0\n",
    "        \n",
    "        else:\n",
    "            vh_array = pad_image(f.ReadAsArray(), tilesize)\n",
    "            invalid_pixels = np.nonzero(vh_array == 0.0)\n",
    "            vh_tiles = tile_image(vh_array,width=tilesize,height=tilesize)\n",
    "            a = np.shape(vh_tiles)\n",
    "            vh_std = np.zeros(a[0])\n",
    "            vht_masked = np.ma.masked_where(vh_tiles<=0, vh_tiles)\n",
    "            vh_picktiles = np.zeros_like(vh_tiles)\n",
    "            for k in range(a[0]):\n",
    "                vh_subtiles = tile_image(vht_masked[k,:,:],width=tilesize2,height=tilesize2)\n",
    "                vhst_mean = np.ma.mean(vh_subtiles, axis=(1,2))\n",
    "                vhst_std = np.ma.std(vhst_mean)\n",
    "                vh_std[k] = np.ma.std(vhst_mean)\n",
    "            \n",
    "            # find tiles with largest standard deviations\n",
    "            vh_mean = np.ma.median(vht_masked, axis=(1,2))\n",
    "            x_vh = np.sort(vh_std/vh_mean)\n",
    "            xm_vh = np.sort(vh_mean)\n",
    "            #x_vh = np.sort(vh_std)            \n",
    "            y_vh = np.arange(1, x_vh.size+1) / x_vh.size\n",
    "            ym_vh = np.arange(1, xm_vh.size+1) / xm_vh.size\n",
    "            \n",
    "            percentile2 = precentile\n",
    "            noconverge_vh = 0\n",
    "            sort_index = 0\n",
    "            while np.size(sort_index) < 5 and noconverge_vh == 0: \n",
    "                threshold_index_vh = np.ma.min(np.where(y_vh>percentile2))\n",
    "                threshold_vh = x_vh[threshold_index_vh]\n",
    "                #sd_select_vh = np.nonzero(vh_std/vh_mean>threshold_vh)\n",
    "                s_select_vh = np.nonzero(vh_std/vh_mean>threshold_vh) \n",
    "                h_select_vh = np.nonzero(Hpercent > Hpick)               # Includes HAND-EM in selection\n",
    "                sd_select_vh = np.intersect1d(s_select_vh, h_select_vh)\n",
    "    \n",
    "                # find tiles with mean values lower than the average mean \n",
    "                omean_vh = np.ma.median(vh_mean[h_select_vh])\n",
    "                mean_select_vh = np.nonzero(vh_mean<omean_vh)\n",
    "            \n",
    "                # Intersect tiles with large std with tiles that have small means\n",
    "                msdselect_vh = np.intersect1d(sd_select_vh, mean_select_vh)\n",
    "                sort_index = np.flipud(np.argsort(vh_std[msdselect_vh]))\n",
    "                percentile2 = percentile2 - 0.01\n",
    "                if percentile2 < 0.5:\n",
    "                    noconverge_vh = 1\n",
    "                    sort_index = threshold_index_vh\n",
    "                    m_thresh_vh = (vh_corr/10.0+30)\n",
    "                    print(\"Tile Selection Did Not Converge - Going with Default Threshold\")\n",
    "            \n",
    "            if noconverge_vh < 1:\n",
    "                finalselect_vh = sort_index[0:5]\n",
    "    \n",
    "            temp = np.ma.masked_where(vh_array<=0, vh_array)\n",
    "            np.seterr(divide='ignore')\n",
    "            np.seterr(invalid='ignore')\n",
    "            dbvh = np.log10(temp)+30\n",
    "            \n",
    "            if noconverge_vh < 1:\n",
    "                # find local thresholds for 5 \"best\" tiles in the image\n",
    "                l_thresh_vh = np.zeros(5)\n",
    "                EMthresh_vh = np.zeros(5)\n",
    "                scaling = 256/(np.mean(dbvh) + 3*np.std(dbvh))\n",
    "                #scaling = 256/(np.mean(vh_array) + 3*np.std(vh_array))\n",
    "                np.seterr(divide='ignore')\n",
    "                np.seterr(invalid='ignore')\n",
    "                dbtile = np.log10(vht_masked)+30\n",
    "                for k in range(5):\n",
    "                    test = dbtile[msdselect_vh[finalselect_vh[k]]] * scaling\n",
    "                    #test = vht_masked[msdselect_vh[finalselect_vh[k]]] * scaling\n",
    "                    A = np.around(test)\n",
    "                    A = A.astype(int)\n",
    "                    #t_thresh = Kittler(A)\n",
    "                    [posterior, cm, cv, cp] = EMSeg_opt(A, 3)\n",
    "                    sorti = np.argsort(cm)\n",
    "                    cms = cm[sorti]\n",
    "                    cvs = cv[sorti]\n",
    "                    cps = cp[sorti]\n",
    "                    xvec = np.arange(cms[0],cms[1],step=.05)\n",
    "                    x1 = make_distribution(cms[0], cvs[0], cps[0], xvec)\n",
    "                    x2 = make_distribution(cms[1], cvs[1], cps[1], xvec)\n",
    "                    dx = np.abs(x1 - x2)\n",
    "                    diff1 = posterior[:,:,0] - posterior[:,:,1]\n",
    "                    t_ind = np.argmin(dx)\n",
    "                    EMthresh_vh[k] = xvec[t_ind]/scaling\n",
    "\n",
    "                    #l_thresh_vh[k] = t_thresh / scaling\n",
    "\n",
    "\n",
    "                    # Mark Tiles used for Threshold Estimation\n",
    "                    vh_picktiles[msdselect_vh[finalselect_vh[k]],:,:]= np.ones_like(vh_tiles[msdselect_vh[finalselect_vh[k]],:,:])\n",
    "    \n",
    "            if noconverge_vh < 1:\n",
    "                # Calculate best threshold for VV and VH as the mean of the 5 thresholds calculated in the previous section \n",
    "                #m_thresh_vh = np.median(l_thresh_vh)\n",
    "                #print(EMthresh_vh-30)\n",
    "                EMts = np.sort(EMthresh_vh)\n",
    "                #m_thresh_vh = np.median(EMthresh_vh)\n",
    "                m_thresh_vh = np.median(EMts[0:4])\n",
    "            \n",
    "            print(\"Best VH Flood Mapping Threshold [dB]: %.2f\" % (10*(m_thresh_vh-30)))\n",
    "            print(\" \")\n",
    "    \n",
    "            # Derive flood mask using the best threshold\n",
    "            maskedarray = np.ma.masked_where(dbvh<=0, dbvh)\n",
    "            \n",
    "            #maskedarray = np.ma.masked_where(vh_array==0, vh_array)\n",
    "            if m_thresh_vh < (vh_corr/10.0+30):\n",
    "                change_mag_mask_vh = maskedarray < m_thresh_vh\n",
    "                vh_thresholds_corr = np.append(vh_thresholds_corr, 10.0*(m_thresh_vh-30)) \n",
    "                #change_mag_mask_vv = np.ma.masked_where(vv_array==0, vv_array) < m_thresh_vv\n",
    "            else:\n",
    "                change_mag_mask_vh = maskedarray < (vh_corr/10.0+30)\n",
    "                vh_thresholds_corr = np.append(vh_thresholds_corr, vh_corr) \n",
    "            # change_mag_mask_vh = vh_array < m_thresh_vh\n",
    "    \n",
    "            # Create Binary masks showing flooded pixels as \"1\"s\n",
    "            sel = np.ones_like(vh_array)\n",
    "            flood_vh = np.zeros_like(vh_array)\n",
    "            flood_vh[change_mag_mask_vh] = sel[change_mag_mask_vh]\n",
    "            np.putmask(flood_vh,vh_array==0 , 0)\n",
    "\n",
    "            # Export flood maps as GeoTIFFs\n",
    "            filename, ext = os.path.basename(tiff).split('.')\n",
    "            outfile = f\"{mask_directory}/{filename}_water_mask.{ext}\"\n",
    "            write_mask_to_file(flood_vh, outfile, f.GetProjection(), f.GetGeoTransform())\n",
    "            vh_array = 0\n",
    "            maskedarray = 0\n",
    "            temp = 0\n",
    "            dbvh = 0\n",
    "            change_mag_mask_vh = 0\n",
    "            Hmask = 0\n",
    "            Hpercent = 0\n",
    "            Hsum = 0\n",
    "            hand_tiles = 0\n",
    "            handem = 0\n",
    "            handem_p = 0\n",
    "            vh_tiles = 0\n",
    "            \n",
    "        \n",
    "    \n",
    "    # Create Maps (Pickfiles) that show which tiles were used for adaptive threshold calculation\n",
    "    if noconverge_vv < 1:\n",
    "        vv_picks = vv_picktiles.reshape((n_rows, n_cols, tilesize, tilesize)) \\\n",
    "                    .swapaxes(1, 2) \\\n",
    "                    .reshape(n_rows * tilesize, n_cols * tilesize)  # yapf: disable\n",
    "    if noconverge_vh < 1:\n",
    "        vh_picks = vh_picktiles.reshape((n_rows, n_cols, tilesize, tilesize)) \\\n",
    "                    .swapaxes(1, 2) \\\n",
    "                    .reshape(n_rows * tilesize, n_cols * tilesize)  # yapf: disable\n",
    "    vh_picktiles = 0\n",
    "    vv_picktiles = 0\n",
    "    \n",
    "    # Write Pickfiles to GeoTIFFs\n",
    "    #if noconverge_vv < 1:\n",
    "    #    outfile = f\"{mask_directory}/{filename[:-3]}_vv_pickfile.{ext}\"\n",
    "    #    write_mask_to_file(vv_picks, outfile, f.GetProjection(), f.GetGeoTransform())\n",
    "    #if noconverge_vh < 1:\n",
    "    #    outfile = f\"{mask_directory}/{filename[:-3]}_vh_pickfile.{ext}\"\n",
    "    #    write_mask_to_file(vh_picks, outfile, f.GetProjection(), f.GetGeoTransform())\n",
    "\n",
    "    # Combine VV and VH flood maps to produce a combined flood mapping product\n",
    "    comb = flood_vh + flood_vv\n",
    "    comb_mask = comb > 0\n",
    "    flood_comb = np.zeros_like(vv_array)\n",
    "    flood_comb[comb_mask] = sel[comb_mask]\n",
    "    filename, ext = os.path.basename(tiff).split('.')\n",
    "    outfile = f\"{mask_directory}/{filename[:-3]}_water_mask_combined.{ext}\"\n",
    "    write_mask_to_file(flood_comb, outfile, f.GetProjection(), f.GetGeoTransform())\n",
    "    \n",
    "    # Create Information on Thresholds used as well as Flood extent information in km2\n",
    "    vv_thresholds = np.append(vv_thresholds, 10.0*(m_thresh_vv-30))\n",
    "    vh_thresholds = np.append(vh_thresholds, 10.0*(m_thresh_vh-30)) \n",
    "    floodarea = np.append(floodarea,(np.sum(flood_comb)*30**2./(1000**2)))\n",
    "    flood_vh = 0\n",
    "    flood_vv = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Post Processing Steps to Clean up Initial Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Use Fuzzy Logic to Improve the Quality of The Water Mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=3>This step applies fuzzy logic roles to remove spurious false detection and improve upon the initial flood mapping product created in Section 2.5. Four Fuzzy indicators for the presence of floods are used. These include:\n",
    "<ol>\n",
    "    <li>The radar cross section in a pixel relative to the determined detection threshold.</li>\n",
    "    <li>The HAND elevation (surface elevation relative to the nearest dranage basin).</li>\n",
    "    <li>The surface slope.</li>\n",
    "    <li>The size of a detected flood feature.</li>\n",
    "</ol>\n",
    "Fuzzy membership functions are calculated for each of these four indicators using a Z-shaped activation function. Membership functions are combined using arithmetric averaging to form a final decision map. This map is then thresholded using a fuzzy threshold of 0.45. <br><br>\n",
    "\n",
    "<b><u>Upper and lower thresholds for the fuzzy activation functions are calculated as follows:</u></b>\n",
    "<ol>\n",
    "    <li><b>RCS:</b> $\\begin{Bmatrix} \n",
    "    x_{u,RCS} & = & \\tau_g \\\\\n",
    "    x_{l,RCS} & = & \\mu_{\\sigma^0(\\tau_g)}\n",
    "    \\end{Bmatrix}$ \n",
    "    with $\\sigma^0(\\tau_g)$ = initial flood classification and flood mapping threshold $\\tau_g$.<br><br></li>\n",
    "    <li><b>HAND:</b> $\\begin{Bmatrix} \n",
    "    x_{u,HAND} & = & \\mu_{HAND(water)} + 3 \\cdot \\sigma_{HAND(water)} \\\\\n",
    "    x_{l,HAND} & = & \\mu_{HAND(water)}\n",
    "    \\end{Bmatrix}$. This is a departure from <cite><a href=\"https://www.sciencedirect.com/science/article/pii/S0924271614001981\" target=\"_blank\"><i>the paper</i></a> by Martinis et al. (2015)</cite>.<br><br></li>\n",
    "    <li><b>Surface slope $\\alpha$:</b> $\\begin{Bmatrix} \n",
    "    x_{u,\\alpha} & = & 0^{\\circ} \\\\\n",
    "    x_{l,\\alpha} & = & 15^{\\circ}\n",
    "    \\end{Bmatrix}$.<br><br></li>\n",
    "    <li><b>Area $A$:</b>$\\begin{Bmatrix} \n",
    "    x_{u,A} & = & 10px \\\\\n",
    "    x_{l,A} & = & 3px\n",
    "    \\end{Bmatrix}$. These threshold values are different differently than in <cite><a href=\"https://www.sciencedirect.com/science/article/pii/S0924271614001981\" target=\"_blank\"><i>the paper</i></a> by Martinis et al. (2015)</cite>.</li>\n",
    "</ol>\n",
    "</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Perform Fuzzy Logic Post Processing of the Initial Flood Maps\n",
    "import skimage.measure\n",
    "import skimage.color\n",
    "from skimage import morphology\n",
    "from skimage import filters\n",
    "import time\n",
    "from IPython.display import clear_output\n",
    "vvcount = 0\n",
    "vhcount = 0\n",
    "\n",
    "print('FUZZY LOGIC POST PROCESSING TO REFINE INITIAL SURFACE WATER EXTENT MAPS:')\n",
    "print(' ')\n",
    "\n",
    "for pair in grouped_pths:\n",
    "    for tiff in grouped_pths[pair]:\n",
    "        start = time.time()\n",
    "        print('-----------------------------------------------------------------------------------')\n",
    "        print(f'Image: {tiff}')\n",
    "        resampled_dem_path=f'{tiff_dir}/resamp_dem.tif'\n",
    "        f = gdal.Open(tiff)\n",
    "        img_array = f.ReadAsArray()\n",
    "        original_shape = img_array.shape\n",
    "        img_array = 0     # free up RAM\n",
    "        radar_array = pad_image(f.ReadAsArray(), tilesize)\n",
    "        temp = np.ma.masked_where(radar_array==0, radar_array)\n",
    "        np.seterr(divide='ignore')\n",
    "        np.seterr(invalid='ignore')\n",
    "        dbarray = np.log10(temp)+30\n",
    "        del temp\n",
    "        \n",
    "        # ------------------------------------------------------------#\n",
    "        # Loading Water Mask File                                     #\n",
    "        #-------------------------------------------------------------#\n",
    "        filename, ext = os.path.basename(tiff).split('.')\n",
    "        waterfile = f\"{mask_directory}/{filename}_water_mask.{ext}\"\n",
    "        h = gdal.Open(waterfile)\n",
    "        maskimage = h.ReadAsArray()\n",
    "    \n",
    "        \n",
    "        print('   - Extracting Relevant Subset from HAND Layer ...')\n",
    "\n",
    "        info1 = (gdal.Info(tiff, options = ['-json']))\n",
    "        info1 = json.dumps(info1)\n",
    "        ul = (json.loads(info1))['cornerCoordinates']['upperLeft']\n",
    "        lr = (json.loads(info1))['cornerCoordinates']['lowerRight']\n",
    "        coordsys = (json.loads(info1))['coordinateSystem']['wkt']\n",
    "        Szone = coordsys.split('ID')[-1].split(',')[1][0:-2]\n",
    "        Sproj_type = coordsys.split('ID')[-1].split('\"')[1]\n",
    "\n",
    "        west= ul[0]\n",
    "        east= lr[0]\n",
    "        south= lr[1]\n",
    "        north= ul[1]\n",
    "        \n",
    "        cmd_resamp=f\"gdalwarp -overwrite -s_srs {Hproj_type}:\"\\\n",
    "        f\"{Hzone} -t_srs EPSG:{Szone} -te {west} {south} {east} {north} -ts {original_shape[1]} {original_shape[0]} -r lanczos {HAND_file} {resampled_dem_path}\"\n",
    "        #print(cmd_resamp)\n",
    "        os.system(cmd_resamp)\n",
    "        \n",
    "        g = gdal.Open(resampled_dem_path)\n",
    "        hand = pad_image(g.ReadAsArray(), tilesize)\n",
    "        \n",
    "        print('   - Interpolate NaNs in HAND Layer ...')\n",
    "        \n",
    "        hand_interp= fill_nan(hand)\n",
    "        del hand\n",
    "        \n",
    "        print('   - Calculate DEM (HAND) Slope magnitude ...')\n",
    "        \n",
    "        vgrad = np.gradient(hand_interp)\n",
    "        mag = np.sqrt(vgrad[0]**2 + vgrad[1]**2)\n",
    "        geotransform = g.GetGeoTransform()\n",
    "        res = geotransform[1]\n",
    "        slope = np.arctan(mag/res)/np.pi*180\n",
    "\n",
    "        print('   - Segment initial flood mask to calculate area of connected patches ...')\n",
    "        \n",
    "        # Here an attempt to perform a sequence of opening and closing steps to \n",
    "        # reduce the number of segments in the inital flood maps and speed up the next processing steps\n",
    "        med = filters.median(maskimage, morphology.disk(2))\n",
    "        selem = morphology.disk(3)\n",
    "        closed = skimage.morphology.closing(med, selem)\n",
    "        #labeled_image = skimage.measure.label(closed, connectivity=2)\n",
    "        labeled_image = skimage.measure.label(maskimage, connectivity=2)\n",
    "        label_areas = np.bincount(labeled_image.ravel())[1:]\n",
    "        \n",
    "        # Define upper and lower limit of the Z fuzzy activation function\n",
    "        \n",
    "        print('   - Now Calculate Fuzzy Weights ...')\n",
    "        \n",
    "        ## sigma zero upper and lower fuzzy threshold calculation\n",
    "        maskedarray = np.ma.masked_where((maskimage==0) | (radar_array < 0), radar_array)\n",
    "        db_lowerlimit = np.log10(np.ma.median(maskedarray))+30\n",
    "        if 'vv' in tiff or 'VV' in tiff:\n",
    "            db_upperlimit = (vv_thresholds_corr[vvcount]/10.0+30)\n",
    "            vvcount = vvcount + 1\n",
    "        else:\n",
    "            db_upperlimit = (vh_thresholds_corr[vhcount]/10.0+30)\n",
    "            vhcount = vhcount + 1\n",
    "        \n",
    "        ## HAND upper and lower fuzzy threshold calculation\n",
    "        maskedarray = np.ma.masked_where(maskimage==0, hand_interp)\n",
    "        ma2 = np.ma.masked_where(maskedarray > np.percentile(maskedarray, 90), maskedarray)\n",
    "        hand_lowerlimit = np.ma.median(np.ma.masked_invalid(ma2))\n",
    "        #hand_upperlimit = hand_lowerlimit + (np.ma.std(np.ma.masked_invalid(ma2)) + 4.5)*np.ma.std(np.ma.masked_invalid(ma2))\n",
    "        hand_upperlimit = hand_lowerlimit + 3.0*(np.ma.std(np.ma.masked_invalid(ma2))) + 5.0\n",
    "\n",
    "        # Create vector spanning all possible values for db, HAND, Slope, and Area values for the selected data set\n",
    "        x_db = np.arange(np.min(dbarray), np.max(dbarray), 0.005)\n",
    "        x_hand = np.arange(np.min(np.ma.masked_invalid(hand_interp)), np.max(np.ma.masked_invalid(hand_interp)), 0.1)\n",
    "        x_slope = np.arange(np.min(np.ma.masked_invalid(slope)), np.max(np.ma.masked_invalid(slope)), 0.1)\n",
    "\n",
    "        largestCC = labeled_image == np.argmax(np.bincount(labeled_image.flat, weights=maskimage.flat))\n",
    "        x_area = np.arange(1, np.sum(largestCC)+10, 1)\n",
    "\n",
    "        # Create activation functions for db, HAND, SLope and Area\n",
    "        activation_db = skfuzzy.zmf(x_db,db_lowerlimit,db_upperlimit)\n",
    "        activation_hand = skfuzzy.zmf(x_hand,hand_lowerlimit,hand_upperlimit)\n",
    "        activation_slope = skfuzzy.zmf(x_slope,0,15)\n",
    "        activation_area = 1-skfuzzy.zmf(x_area,3,10)\n",
    "        \n",
    "        # Calculate membership functions given your activation function rule\n",
    "        db_membership = skfuzzy.interp_membership(x_db, activation_db, dbarray)\n",
    "        print('      -- Calculating Radar Cross Section Membership ... [min thresholds: %6.2f' % (db_lowerlimit-30),'; max threshold: %6.2f]' % (db_upperlimit-30))     \n",
    "        hand_membership = skfuzzy.interp_membership(x_hand, activation_hand, hand_interp)\n",
    "        print('      -- Calculating HAND Membership ... [min thresholds: %6.2f' % (hand_lowerlimit),'; max threshold: %6.2f]' % (hand_upperlimit)) \n",
    "        slope_membership = skfuzzy.interp_membership(x_slope, activation_slope, slope)\n",
    "        print('      -- Calculating Slope Membership ... [min thresholds: 0 deg; max threshold: 15 deg]') \n",
    "        area_membership = np.zeros_like(radar_array)\n",
    "        \n",
    "        # ALM\n",
    "        # try to speed up 'area_membership'\n",
    "        # first use np.histogram to create simple arrays of label numbers and corresponding areas (counts of pixels)\n",
    "        # in this manner the 'bin' center is the label number and the 'count' is the area for that label\n",
    "        area_membership = np.zeros_like(radar_array)\n",
    "        nlabels = np.max(labeled_image)\n",
    "        bin_edges = np.arange(nlabels+2)-0.5\n",
    "        label_sizes, bin_edges = np.histogram(labeled_image, bins=bin_edges)\n",
    "        \n",
    "        # ALM\n",
    "        # the fuzzy membership for the activation area is 0 for areas less than area_min (3)\n",
    "        # the fuzzy membership for the activation area is 1 for areas greater than area_max (10)\n",
    "        # therefore we only really need to calculate activation area for a discrete set of values\n",
    "        area_min = 3\n",
    "        area_max = 10\n",
    "        \n",
    "        # find areas below the min\n",
    "        import time\n",
    "        start_time = time.time()\n",
    "        print('starting area membership -- ALM:')\n",
    "        match = np.where(label_sizes < area_min)\n",
    "        mask = np.isin(labeled_image, match)\n",
    "        np.putmask(area_membership, mask, 0)\n",
    "        # find areas above the max\n",
    "        match = np.where(label_sizes > area_max)\n",
    "        mask = np.isin(labeled_image, match)\n",
    "        np.putmask(area_membership, mask, 1)\n",
    "        # iterate through the rest of them\n",
    "        areas = np.arange(area_min, area_max+1)\n",
    "        for area in areas:\n",
    "            match = np.where(label_sizes == area)\n",
    "            mask = np.isin(labeled_image, match)\n",
    "            wt = skfuzzy.interp_membership(x_area, activation_area, area)\n",
    "            np.putmask(area_membership, mask, wt)\n",
    "        print('finished area membership -- ALM')\n",
    "        print('--- %s seconds ---' % (time.time()-start_time))\n",
    "#       ALM new contributions end here ^^^\n",
    "        \n",
    "        # print('      -- Calculating Area Membership ... [This may take a while!]')\n",
    "        # for x in tqdm(range(1, np.max(labeled_image))):\n",
    "            #clear_output(wait=True)\n",
    "            #np.putmask(area_membership,labeled_image==x,skfuzzy.interp_membership(x_area, activation_area, np.sum(labeled_image==x)))\n",
    "            # np.putmask(area_membership,labeled_image==x,skfuzzy.interp_membership(x_area, activation_area, label_areas[x-1]))\n",
    "            \n",
    "        water_gT=gdal_get_geotransform(waterfile)\n",
    "        water_proj4=get_proj4(waterfile)\n",
    "        \n",
    "        # Uncomment these following lines if you want more intermediary files to be saved off as GeoTIFFs\n",
    "            #filename, ext = os.path.basename(tiff).split('.')\n",
    "            #outfile = f\"{mask_directory}/{filename}_dbmembership.{ext}\"\n",
    "            #gdal_write(db_membership, water_gT, filename=outfile, srs_proj4=water_proj4.srs, nodata=np.nan, format=gdal.GDT_Float32)\n",
    "            #outfile = f\"{mask_directory}/{filename}_handmembership.{ext}\"\n",
    "            #gdal_write(hand_membership, water_gT, filename=outfile, srs_proj4=water_proj4.srs, nodata=np.nan, format=gdal.GDT_Float32)\n",
    "            #outfile = f\"{mask_directory}/{filename}_slopemembership.{ext}\"\n",
    "            #gdal_write(slope_membership, water_gT, filename=outfile, srs_proj4=water_proj4.srs, nodata=np.nan, format=gdal.GDT_Float32)\n",
    "            #outfile = f\"{mask_directory}/{filename}_areamembership.{ext}\"\n",
    "            #gdal_write(area_membership, water_gT, filename=outfile, srs_proj4=water_proj4.srs, nodata=np.nan, format=gdal.GDT_Float32)\n",
    "        \n",
    "       \n",
    "        print('   - Combine all fuzzy membership functions to make final flood mapping decision')\n",
    "        dbmask = db_membership != 0.0\n",
    "        handmask = hand_membership != 0.0\n",
    "        slopemask = slope_membership != 0.0\n",
    "        areamask = area_membership != 0.0\n",
    "        combinedm = dbmask * handmask * slopemask * areamask\n",
    "        combinedmask = np.zeros_like(radar_array)\n",
    "        sel = np.ones_like(dbarray)\n",
    "        combinedmask[combinedm] = sel[combinedm]\n",
    "        combinedweights = ((db_membership + hand_membership + slope_membership + area_membership) / 4.0) * combinedmask \n",
    "        acceptance = combinedweights > 0.45\n",
    "        \n",
    "        # Uncomment these following lines if you want more intermediary files to be saved off as GeoTIFFs\n",
    "        #outfile = f\"{mask_directory}/{filename}_totalmembership.{ext}\"\n",
    "        #gdal_write(combinedweights, water_gT, filename=outfile, srs_proj4=water_proj4.srs, nodata=np.nan, format=gdal.GDT_Float32)\n",
    "        \n",
    "        sel = np.ones_like(radar_array)\n",
    "        floodmap = np.zeros_like(radar_array)\n",
    "        floodmap[acceptance] = sel[acceptance]\n",
    "        np.putmask(floodmap,radar_array==0 , 0)\n",
    "        if 'vv' in tiff or 'VV' in tiff:\n",
    "            floodmap_vv = floodmap\n",
    "        else:\n",
    "            floodmap_vh = floodmap\n",
    "\n",
    "        # Export flood maps as GeoTIFFs\n",
    "        #filename, ext = os.path.basename(tiff).split('.')\n",
    "        #outfile = f\"{mask_directory}/{filename}_final_water_mask.{ext}\"\n",
    "        #write_mask_to_file(floodmap, outfile, f.GetProjection(), f.GetGeoTransform())\n",
    "        del floodmap\n",
    "        print('   - Processing time: %6.2f minutes' % ((time.time() - start)/60.0))\n",
    "        print('-----------------------------------------------------------------------------------')\n",
    "        print(' ')\n",
    "    \n",
    "    # Combine VV and VH flood maps to produce a combined flood mapping product\n",
    "    comb = floodmap_vh + floodmap_vv\n",
    "    comb_mask = comb > 0\n",
    "    flood_comb = np.zeros_like(radar_array)\n",
    "    flood_comb[comb_mask] = sel[comb_mask]\n",
    "    filename, ext = os.path.basename(tiff).split('.')\n",
    "    outfile = f\"{mask_directory}/{filename[:-3]}_fcWM.{ext}\"\n",
    "    write_mask_to_file(flood_comb, outfile, f.GetProjection(), f.GetGeoTransform())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Version Log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"2\" color=\"gray\"> <i> Version 1.01 - Franz J Meyer; 12/15/2020\n",
    "\n",
    "Recent Changes:\n",
    "- Post processing was integrated\n",
    "- Notebook was trimmed down to the necessary parts\n",
    "- Analysis of results moved into a separate notebook\n",
    "</i></font>"
   ]
  }
 ],
 "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.8.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}