Flood_Depth_Estimation_with_Flood_Extent_Maps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<img src=\"https://radar.community.uaf.edu/wp-content/uploads/sites/667/2021/03/HydroSARbanner.jpg\" width=\"100%\" />\n",
    "<br>\n",
    "<font face=\"Calibri\">\n",
    "\n",
    "<font size=\"6\"> <b> Flood Depth Estimation with Flood Extent Maps </b> </font> \n",
    "\n",
    "<font size=\"4\"> <b> Part of NASA A.37 Project:</b> Integrating SAR Data for Improved Resilience and Response to Weather-Related Disasters   <br>\n",
    "<font size=\"4\"> <b> PI:</b>Franz J. Meyer <br>\n",
    "<font size=\"3\"> Version 0.1.8 - 2021/01/24 <br>\n",
    "<b>Change Log:</b> See bottom of the notebook.<br>\n",
    "    \n",
    "<br>\n",
    "<font size=\"4\"> <b> Batuhan Osmanoglu, MinJeong Jo; NASA Goddard Space Fligth Center </b> <br>\n",
    "</font>\n",
    "    \n",
    "<font size=\"3\"> This notebook provides the processor to generate Flood Depth map using the product generated by <b><font color='rgba(200,0,0,0.2)'>Hyp3 Change Detection-Threshold </font></b> processor. This notebook can be used to generate <b><font color='rgba(200,0,0,0.2)'> Multiple </font></b> FD Products  </font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\">\n",
    "<font size=\"3\"><b>Note:</b> Before you start to use the notebook, <b><font color='rgba(200,0,0,0.2)'>Hyp3-generated change detection maps in Geotiff format </font></b> need to be placed in your own data folder. <b><font color='rgba(200,0,0,0.2)'>For the HydroSAR Training, these SAR data are already available to you after completion of Lab 2.  </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "# Importing Relevant Python Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#Setup Environment\n",
    "import os\n",
    "import urllib\n",
    "import numpy as np\n",
    "import gdal\n",
    "import osr\n",
    "import json\n",
    "import glob\n",
    "import scipy # in first try, had scipy error\n",
    "from scipy import ndimage\n",
    "from scipy import interpolate\n",
    "from scipy import optimize\n",
    "from scipy import stats\n",
    "from importlib import reload\n",
    "from skimage import morphology\n",
    "import pylab as pl\n",
    "import pprint\n",
    "import concurrent.futures as cf\n",
    "from IPython.core.debugger import set_trace\n",
    "#The two lines below are for visually browsing and selecting the DEM. \n",
    "import ipywidgets as ui\n",
    "from IPython.display import display\n",
    "\n",
    "import warnings\n",
    "\n",
    "# try:\n",
    "#     import cv2\n",
    "# except:\n",
    "#     !pip install --user opencv-python\n",
    "#    import cv2\n",
    "try:\n",
    "    import astropy\n",
    "except:\n",
    "    !pip install --user astropy\n",
    "    import astropy\n",
    "import astropy.convolution\n",
    "try:\n",
    "    import pykrige\n",
    "except:\n",
    "    !pip install --user pykrige\n",
    "    import pykrige\n",
    "try:\n",
    "    import pysheds\n",
    "except:\n",
    "    !pip install pysheds\n",
    "    import pysheds\n",
    "from pysheds.grid import Grid\n",
    "try:\n",
    "    from affine import Affine\n",
    "except:\n",
    "    !pip install affine\n",
    "    from affine import Affine\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",
    "#Download packages\n",
    "codes_folder='/home/jovyan/codes'\n",
    "project_dir=os.getcwd()\n",
    "try:\n",
    "    os.chdir(codes_folder)\n",
    "except:\n",
    "    os.mkdir(codes_folder)\n",
    "    os.chdir(codes_folder)\n",
    "\n",
    "if not os.path.exists(codes_folder+'/adore-doris/'):\n",
    "  !git clone https://github.com/bosmanoglu/adore-doris.git\n",
    "os.sys.path.append(codes_folder+'/adore-doris/lib/python')\n",
    "os.sys.path.append(codes_folder)\n",
    "os.chdir(project_dir) #Go back to project folder\n",
    "\n",
    "#import modules after downloads\n",
    "import gis\n",
    "import surface_from_grad as sfg\n",
    "import deramp as surface\n",
    "import invdisttree\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "#%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Define Convenience Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     1,
     17,
     21,
     29,
     33,
     37,
     62,
     70,
     73,
     89,
     98,
     113,
     119,
     125,
     141,
     158,
     194,
     293,
     310,
     323,
     336
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "# Define convenience functions\n",
    "def bounding_box_inside_bounding_box(small, big):\n",
    "    s0=np.array([p[0] for p in small])\n",
    "    s1=np.array([p[1] for p in small])\n",
    "    b0=np.array([p[0] for p in big])\n",
    "    b1=np.array([p[1] for p in big])\n",
    "    inside=True\n",
    "    if s0.min()<b0.min():\n",
    "        inside=False\n",
    "    if s0.max()>b0.max():\n",
    "        inside=False\n",
    "    if s1.min()<b1.min():\n",
    "        inside=False\n",
    "    if s1.max()>b1.max():\n",
    "        inside=False\n",
    "    return inside\n",
    "\n",
    "def getGeoTransform(filename):\n",
    "    warnings.warn(\"getGeoTransform will be deprecated in the future. Please use read_data instead.\", PendingDeprecationWarning)\n",
    "    return get_geotransform(filename)\n",
    "\n",
    "def get_geotransform(filename):\n",
    "    '''\n",
    "    [top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution]=getGeoTransform('/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 build_vrt(filename, input_file_list):\n",
    "    vrt_options = gdal.BuildVRTOptions(resampleAlg='near', resolution='highest', separate=False, targetAlignedPixels=True)\n",
    "    gdal.BuildVRT(filename,input_file_list,options=vrt_options)\n",
    "\n",
    "def get_tiff_paths(paths):\n",
    "    tiff_paths = !ls $paths | sort -t_ -k5,5\n",
    "    return tiff_paths\n",
    "\n",
    "def gdal_get_projection(filename, out_format='proj4'):\n",
    "    \"\"\"\n",
    "    epsg_string=get_epsg(filename, out_format='proj4')\n",
    "    \"\"\"\n",
    "    try:\n",
    "      ds=gdal.Open(filename, gdal.GA_ReadOnly)\n",
    "      srs=gdal.osr.SpatialReference()\n",
    "      srs.ImportFromWkt(ds.GetProjectionRef())\n",
    "    except: #I am not sure if this is working for datasets without a layer. The first try block should work mostly.\n",
    "      ds=gdal.Open(filename, gdal.GA_ReadOnly)\n",
    "      ly=ds.GetLayer()\n",
    "      if ly is None:\n",
    "        print(f\"Can not read projection from file:{filename}\")\n",
    "        return None\n",
    "      else:\n",
    "        srs=ly.GetSpatialRef()\n",
    "    if out_format.lower()=='proj4':\n",
    "      return srs.ExportToProj4()\n",
    "    elif out_format.lower()=='wkt':\n",
    "      return srs.ExportToWkt()\n",
    "    elif out_format.lower()=='epsg':\n",
    "      crs=pyproj.crs.CRS.from_proj4(srs.ExportToProj4())\n",
    "      return crs.to_epsg()\n",
    "\n",
    "\n",
    "def get_size(filename):\n",
    "    \"\"\"(width, height) = get_size(filename)\n",
    "    \"\"\"\n",
    "    ds = gdal.Open(filename)\n",
    "    width = ds.RasterXSize\n",
    "    height = ds.RasterYSize\n",
    "    ds=None\n",
    "    return (width, height)\n",
    "def get_proj4(filename):\n",
    "    f=rasterio.open(filename)\n",
    "    return pyproj.Proj(f.crs, preserve_units=True)  #used in pysheds\n",
    "def clip_gT(gT, xmin, xmax, ymin, ymax, method='image'):\n",
    "    '''calculate new geotransform for a clipped raster either using pixels or projected coordinates.\n",
    "    clipped_gT=clip_gT(gT, xmin, xmax, ymin, ymax, method='image')\n",
    "    method: 'image' | 'coord'\n",
    "    '''\n",
    "    if method == 'image':\n",
    "      y,x=xy2coord(ymin, xmin, gT); #top left, reference, coordinate\n",
    "    if method == 'coord':\n",
    "      #find nearest pixel\n",
    "      yi, xi = coord2xy(ymin, xmin, gT)\n",
    "      #get pixel coordinate\n",
    "      y,x=xy2coord(yi, xi, gT)\n",
    "    gTc=list(gT)\n",
    "    gTc[0]=y\n",
    "    gTc[3]=x\n",
    "    return tuple(gTc)\n",
    "def xy2coord(x,y,gT):\n",
    "    '''\n",
    "    lon,lat=xy2coord(x,y,geoTransform)\n",
    "    projects pixel index to position based on geotransform.\n",
    "    '''\n",
    "    coord_x=gT[0] + x*gT[1] + y*gT[2]\n",
    "    coord_y=gT[3] + x*gT[4] + y*gT[5]\n",
    "    return coord_x, coord_y\n",
    "\n",
    "def coord2xy(x,y,gT):\n",
    "    '''\n",
    "    x,y = coord2xy(lon, lat, geoTransform)\n",
    "    calculates pixel index closest to the lon, lat.\n",
    "    '''\n",
    "    #ref: https://gis.stackexchange.com/questions/221292/retrieve-pixel-value-with-geographic-coordinate-as-input-with-gdal/221430\n",
    "    xOrigin = gT[0]\n",
    "    yOrigin = gT[3]\n",
    "    pixelWidth = gT[1]\n",
    "    pixelHeight = -gT[5]\n",
    "\n",
    "    col = np.array((x - xOrigin) / pixelWidth).astype(int)\n",
    "    row = np.array((yOrigin - y) / pixelHeight).astype(int)\n",
    "\n",
    "    return row,col\n",
    "def fitSurface(x,y,z,X,Y):\n",
    "    p0=[0,0.1,0.1]\n",
    "    fitfunc = lambda p, x, y: p[0]+p[1]*x+p[2]*y\n",
    "    errfunc = lambda p, x, y, z: abs(fitfunc(p,x,y) - z)\n",
    "    planefit, success=optimize.leastsq(errfunc, p0, args=(x,y,z))\n",
    "    return fitfunc(planefit, X,Y)    \n",
    "def nonan(A, rows=False):\n",
    "    if rows:\n",
    "        return A[np.isnan(A).sum(1)==0];\n",
    "    else:\n",
    "        return A[~np.isnan(A)];\n",
    "\n",
    "def get_wesn(filename, t_srs=None):\n",
    "    bb=bounding_box(filename, t_srs=t_srs)\n",
    "    w=np.inf\n",
    "    e=-np.inf\n",
    "    n=-np.inf\n",
    "    s=np.inf\n",
    "    for p in bb:\n",
    "        if p[0]<w:\n",
    "            w=p[0]\n",
    "        if p[0]>e:\n",
    "            e=p[0]\n",
    "        if p[1]<s:\n",
    "            s=p[1]\n",
    "        if p[1]>n:\n",
    "            n=p[1]\n",
    "    return [w,e,s,n]\n",
    "def bounding_box(filename, t_srs=None):\n",
    "    \"\"\"\n",
    "    ((lon1,lat1), (lon2,lat2), (lon3,lat3), (lon4,lat4))=bounding_box('/path/to/file', t_srs=None) #returns x,y in native coordinate system\n",
    "    ((lon1,lat1), (lon2,lat2), (lon3,lat3), (lon4,lat4))=bounding_box('/path/to/file', t_srs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')\n",
    "    \"\"\"\n",
    "    gT=getGeoTransform(filename)\n",
    "    width, height=get_size(filename)\n",
    "    pts=(xy2coord(0,0,gT), xy2coord(width,0,gT), xy2coord(width, height,gT), xy2coord(0, height,gT))\n",
    "    if t_srs is None:\n",
    "        return pts\n",
    "    else:\n",
    "        pts_tsrs=[]\n",
    "        s_srs=gdal_get_projection(filename, out_format='proj4')\n",
    "        for p in pts:\n",
    "            pts_tsrs.append(transform_point(p[0], p[1], 0, s_srs=s_srs, t_srs=t_srs))\n",
    "    return tuple(pts_tsrs)   \n",
    "\n",
    "def transform_point(x,y,z,s_srs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs', t_srs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'):\n",
    "    '''\n",
    "    transform_point(x,y,z,s_srs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs', t_srs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')\n",
    "    \n",
    "    Known Bugs: gdal transform may fail if a proj4 string can not be found for the EPSG or WKT formats. \n",
    "    '''    \n",
    "    srs_cs=osr.SpatialReference()    \n",
    "    if \"EPSG\" == s_srs[0:4]:    \n",
    "      srs_cs.ImportFromEPSG(int(s_srs.split(':')[1]));\n",
    "    elif \"GEOCCS\" == s_srs[0:6]:\n",
    "      srs_cs.ImportFromWkt(s_srs);\n",
    "    else:\n",
    "      srs_cs.ImportFromProj4(s_srs);\n",
    "\n",
    "    trs_cs=osr.SpatialReference()    \n",
    "    if \"EPSG\" == t_srs[0:4]:    \n",
    "      trs_cs.ImportFromEPSG(int(t_srs.split(':')[1]));\n",
    "    elif \"GEOCCS\" == t_srs[0:6]:\n",
    "      trs_cs.ImportFromWkt(t_srs);\n",
    "    else:\n",
    "      trs_cs.ImportFromProj4(t_srs);\n",
    "    if int(gdal.VersionInfo())>2999999: #3010300\n",
    "        #https://gdal.org/tutorials/osr_api_tut.html#crs-and-axis-order\n",
    "        # https://github.com/OSGeo/gdal/blob/master/gdal/MIGRATION_GUIDE.TXT\n",
    "        srs_cs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)\n",
    "        trs_cs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)\n",
    "    transform = osr.CoordinateTransformation(srs_cs,trs_cs) \n",
    "\n",
    "    if numel(x)>1:\n",
    "      return [  transformPoint(x[k], y[k], z[k]) for k in range(numel(x))]\n",
    "    else:\n",
    "      try:\n",
    "        return transform.TransformPoint((x,y,z));\n",
    "      except: \n",
    "        return transform.TransformPoint(x,y,z)\n",
    "\n",
    "def get_waterbody(filename, ths):\n",
    "\n",
    "    corners=bounding_box(filename)\n",
    "\n",
    "    epsg=gdal_get_projection(filename, out_format='epsg')\n",
    "    if epsg==\"4326\":         \n",
    "        corners=bounding_box(filename)\n",
    "    else:\n",
    "        srs=gdal_get_projection(filename, out_format='proj4')\n",
    "        corners=bounding_box(filename, t_srs=\"EPSG:4326\")\n",
    "    west= corners[0][0]\n",
    "    east= corners[1][0]\n",
    "    south= corners[2][1]\n",
    "    north= corners[0][1]    \n",
    "        \n",
    "\n",
    "    cwd = os.getcwd()\n",
    "    sw_path = f\"{cwd}/S_WATER\"\n",
    "    if not os.path.exists(sw_path):\n",
    "        os.mkdir(sw_path)\n",
    "\n",
    "    lon = np.floor(west/10)\n",
    "    lon=int(abs(lon*10))\n",
    "    lat = np.ceil(north/10)\n",
    "    lat=int(abs(lat*10))\n",
    "\n",
    "    if (west<0 and north<0):\n",
    "        urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}W_{lat}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}W_{lat}S.tif\")\n",
    "        if (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon-10}W_{lat}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon-10}W_{lat}S.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}W_{lat+10}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}W_{lat+10}S.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)) and (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon-10}W_{lat+10}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon-10}W_{lat+10}S.tif\")\n",
    "        print(f\"lon: {lon}-{lon-10}W, lat: {lat}-{lat+10}S \")\n",
    "\n",
    "    elif (west<0 and north>=0):\n",
    "        urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}W_{lat}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}W_{lat}N.tif\")\n",
    "        if (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon-10}W_{lat}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon-10}W_{lat}N.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}W_{lat-10}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}W_{lat-10}N.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)) and (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon-10}W_{lat-10}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon-10}W_{lat-10}N.tif\")\n",
    "        print(f\"lon: {lon}-{lon-10}W, lat: {lat}-{lat-10}N \")\n",
    "\n",
    "\n",
    "    elif (west>=0 and north<0):\n",
    "        urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}E_{lat}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}E_{lat}S.tif\")\n",
    "        if (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon+10}E_{lat}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon+10}E_{lat}S.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}E_{lat+10}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}E_{lat+10}S.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)) and (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon+10}E_{lat+10}Sv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon+10}E_{lat+10}S.tif\")\n",
    "        print(f\"lon: {lon}-{lon+10}E, lat: {lat}-{lat+10}S \")\n",
    "\n",
    "    else:\n",
    "        urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}E_{lat}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}E_{lat}N.tif\")\n",
    "        if (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon+10}E_{lat}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon+10}E_{lat}N.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon}E_{lat-10}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon}E_{lat-10}N.tif\")\n",
    "        if (np.floor(north/10) != np.floor(south/10)) and (np.floor(west/10) != np.floor(east/10)):\n",
    "            urllib.request.urlretrieve(f\"https://storage.googleapis.com/global-surface-water/downloads2019v2/occurrence/occurrence_{lon+10}E_{lat-10}Nv1_1_2019.tif\", f\"{cwd}/S_WATER/surface_water_{lon+10}E_{lat-10}N.tif\")\n",
    "        print(f\"lon: {lon}-{lon+10}E, lat: {lat}-{lat-10}N \")\n",
    "\n",
    "    # Building the virtual raster for Change Detection product(tiff)\n",
    "    product_wpath = f\"{cwd}/S_WATER/surface_water*.tif\"\n",
    "    #wildcard_path = f\"{cwd}/change_VV_20170818T122205_20170830T122203.tif\"\n",
    "    print(product_wpath)\n",
    "    os.chdir(cwd)\n",
    "\n",
    "    get_ipython().system('gdalbuildvrt surface_water_map.vrt $product_wpath')\n",
    "\n",
    "    #Clipping/Resampling Surface Water Map for AOI\n",
    "    dim =get_size(filename)\n",
    "    if epsg==\"4326\":\n",
    "        cmd_resamp=f\"gdalwarp -overwrite -te {west} {south} {east} {north} -ts {dim[0]} {dim[1]} -r nearest {cwd}/surface_water_map.vrt {cwd}/surface_water_map_clip.tif\"\n",
    "    else:   \n",
    "        corners=bounding_box(filename) # we now need corners in the non EPSG:4326 format.  \n",
    "        west= corners[0][0]\n",
    "        east= corners[1][0]\n",
    "        south= corners[2][1]\n",
    "        north= corners[0][1]            \n",
    "        cmd_resamp=f\"gdalwarp -overwrite -t_srs '{srs}' -te {west} {south} {east} {north} -ts {dim[0]} {dim[1]} -r nearest {cwd}/surface_water_map.vrt {cwd}/surface_water_map_clip.tif\"        \n",
    "    print(cmd_resamp)\n",
    "    os.system(cmd_resamp)\n",
    "\n",
    "    #load resampled water map\n",
    "    wimage_file =f\"{cwd}/surface_water_map_clip.tif\"\n",
    "    #water_map = gdal.Open(wimage_file)\n",
    "    #print(f\"X-dimension: {water_map.RasterXSize} Y-dimension: {water_map.RasterYSize}\")\n",
    "    \n",
    "    swater_map =gis.readData(wimage_file)\n",
    "    wmask=swater_map>ths   #higher than 30% possibility (present water)\n",
    "        \n",
    "    return wmask\n",
    "\n",
    "def numel(x):\n",
    "    if isinstance(x, np.int):\n",
    "      return 1\n",
    "    elif isinstance(x, np.double):\n",
    "      return 1\n",
    "    elif isinstance(x, np.float):\n",
    "      return 1\n",
    "    elif isinstance(x, str):\n",
    "      return 1\n",
    "    elif isinstance(x, list) or isinstance(x, tuple):\n",
    "      return len(x)\n",
    "    elif isinstance(x, np.ndarray):\n",
    "      return x.size\n",
    "    else: \n",
    "      print('Unknown type {}.'.format(type(x)))\n",
    "      return None\n",
    "\n",
    "def yesno(yes_no_question=\"[y/n]\"):\n",
    "    while True:\n",
    "        # raw_input returns the empty string for \"enter\"\n",
    "        yes = {'yes','y', 'ye'}\n",
    "        no = {'no','n'}\n",
    "\n",
    "        choice = input(yes_no_question+\"[y/n]\").lower()\n",
    "        if choice in yes:\n",
    "            return True\n",
    "        elif choice in no:\n",
    "            return False\n",
    "        else:\n",
    "            print(\"Please respond with 'yes' or 'no'\")\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 logstat(data, func=np.nanstd):\n",
    "    \"\"\" stat=logstat(data, func=np.nanstd)\n",
    "       calculates the statistic after taking the log and returns the statistic in linear scale.\n",
    "       INF values inside the data array are set to nan. \n",
    "       The func has to be able to handle nan values. \n",
    "    \"\"\"\n",
    "    ld=np.log(data)\n",
    "    ld[np.isinf(ld)]=np.nan #requires func to handle nan-data.\n",
    "    st=func(ld)\n",
    "    return np.exp(st)\n",
    "\n",
    "def iterative(hand, extent, water_levels=range(15)):\n",
    "    #accuracy=np.zeros(len(water_levels))    \n",
    "    #for k,w in enumerate(water_levels):\n",
    "    #    iterative_flood_extent=hand<w\n",
    "    #    TP=np.nansum(np.logical_and(iterative_flood_extent==1, extent==1)) #true positive\n",
    "    #    TN=np.nansum(np.logical_and(iterative_flood_extent==0, extent==0)) # True negative\n",
    "    #    FP=np.nansum(np.logical_and(iterative_flood_extent==1, extent==0)) # False positive\n",
    "    #    FN=np.nansum(np.logical_and(iterative_flood_extent==0, extent==1)) # False negative\n",
    "    #    #accuracy[k]=(TP+TN)/(TP+TN+FP+FN) #accuracy \n",
    "    #    accuracy[k]=TP/(TP+FP+FN) #threat score\n",
    "    #best_water_level=water_levels[np.argmax(accuracy)]\n",
    "    \n",
    "    def _goal_ts(w):\n",
    "        iterative_flood_extent=hand<w # w=water level\n",
    "        TP=np.nansum(np.logical_and(iterative_flood_extent==1, extent==1)) #true positive\n",
    "        TN=np.nansum(np.logical_and(iterative_flood_extent==0, extent==0)) # True negative\n",
    "        FP=np.nansum(np.logical_and(iterative_flood_extent==1, extent==0)) # False positive\n",
    "        FN=np.nansum(np.logical_and(iterative_flood_extent==0, extent==1)) # False negative\n",
    "        return 1-TP/(TP+FP+FN) #threat score #we will minimize goal func, hence 1-threat_score.\n",
    "    #bounds=(min(water_levels), max(water_levels))\n",
    "    #opt_res=optimize.minimize(_goal_ts, max(bounds),method='TNC',bounds=[bounds],options={'xtol':0.1, 'scale':1})\n",
    "\n",
    "    class MyBounds(object):\n",
    "        def __init__(self, xmax=[max(water_levels)], xmin=[min(water_levels)] ):\n",
    "            self.xmax = np.array(xmax)\n",
    "            self.xmin = np.array(xmin)\n",
    "        def __call__(self, **kwargs):\n",
    "            x = kwargs[\"x_new\"]\n",
    "            tmax = bool(np.all(x <= self.xmax))\n",
    "            tmin = bool(np.all(x >= self.xmin))\n",
    "            return tmax and tmin\n",
    "    mybounds = MyBounds()\n",
    "    x0=[np.mean(water_levels)]\n",
    "    opt_res=optimize.basinhopping(_goal_ts, x0,niter=10000, niter_success=100, accept_test=mybounds)\n",
    "    if opt_res.message[0]=='success condition satisfied' or opt_res.message[0]=='requested number of basinhopping iterations completed successfully':\n",
    "        best_water_level=opt_res.x[0]\n",
    "    else:        \n",
    "        best_water_level=np.inf # set to inf to mark unstable solution.    \n",
    "    return best_water_level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "class PathSelector():\n",
    "    \"\"\"\n",
    "    Displays a file selection tree. Any file can be selected. \n",
    "    Selected path can be obtained by: PathSelector.accord.get_title(0)\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 = ()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": [],
    "heading_collapsed": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "# Define Some Common Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\">This section allows you to customize how flood depth estimation is performed. <b>The main paramters that users might want to change are</b>:\n",
    "    \n",
    "<ul>\n",
    " <li>\n",
    "<b>Input File Naming Scheme:</b> This is only relevant if you are interested in mosaicking large areas. This gives you the option of either picking initial flood mapping information created in Lab 2 (naming scheme <i>*_water_mask_combined.tiff</i>) or final post-processed flood mapping information (naming scheme <i>*_fcWM.tiff</i>) for flood depth calculation [<font color='rgba(200,0,0,0.2)'><b> for the HydroSAR training, please do not change this variable from its default</b></font>].\n",
    "    </li><br>\n",
    " \n",
    "  <li><b>Estimator:</b> Three different estimation approaches were implemented and are currently being tested by the HydroSAR team: \n",
    "    <ul style=\"list-style-type: lower-alpha; padding-bottom: 0;\">\n",
    "    <li style=\"margin-left:2em\"><b>Iterative:</b> Basin hopping optimization method to match flooded areas to flood depth estimates given the HAND layer. From our current experience, this is the most accurate, but also the most time consuming approach.</li>\n",
    "    <li style=\"margin-left:2em; padding-bottom: 0;\"><b>Normalized Median Absolute Deviation (nmad):</b> Uses a median operator to estimate the variation to increase robustness in the presence of outliers. [<font color='rgba(200,0,0,0.2)'><b> We will use this approach for the HydroSAR training</b></font>].</li>\n",
    "    <li style=\"margin-left:2em; padding-bottom: 0;\"><b>Logstat:</b> This approach calculates mean and standard deviation of HAND heights in the logarithmic domain to improve robustness for very non-Gaussian data distributions.</li>\n",
    "        <li style=\"margin-left:2em; padding-bottom: 0;\"><b>Numpy:</b> Calculates statistics needed in the approach in linear scale. This approach is least robust to outliers and non-Gaussian distributions.</li>\n",
    "    </ul>\n",
    "\n",
    "</ul>\n",
    "</font>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#parameters setup\n",
    "version=\"0.1.8\"\n",
    "water_classes = [1,2,3,4,5] # 1 has to be a water class, 0 is no water Others are optional.\n",
    "pattern=\"*_water_mask_combined.tiff\" #\"filter_*_amp_Classified.tif\"\n",
    "show_plots=True #turn this off for debugging with IPDB\n",
    "water_level_sigma=3 #use 3*std to estimate max. water height (water level) for each object. Used for numpy, nmad,logstat\n",
    "estimator = \"iterative\" # iterative, numpy, nmad or logstat\n",
    "iterative_bounds=[0,15] #only used for iterative\n",
    "output_prefix='' # Output file is created in the same folder as flood extent. A prefix can be added to the filename.\n",
    "known_water_threshold=30 #Threshold for extracting the known water area in percent. \n",
    "if show_plots:\n",
    "    %matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Prepare Data Set for HAND Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<br>\n",
    "<font face=\"Calibri\" size=\"4\"><b>Enter the path to the directory holding your tiffs:</b>. Here we ask you if you want to calculate flood depth across a mosaic or for a single file. [<font color='rgba(200,0,0,0.2)'><b> Please select \"single file\" for the HydroSAR Training Exercise</b></font>]. </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     1
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Here we ask if we are processing a spatial mosaic or a single file. In case of mosaicking multile files, the final water mask's name will be the same with the work directory name \n",
    "if yesno(\"Would you like to mosaic multiple files (e.g. large area coverage from multiple scenes)?\"):\n",
    "    print(f\"Choose one of the water extent files inside the folder. All files with matching the following will be processed: {pattern}\")\n",
    "    file_folder_func=lambda x: os.path.dirname(x)\n",
    "    single_file=False\n",
    "else:    \n",
    "    print(\"Choose your GDAL compatible Classified water extent file using the file browser below:\")\n",
    "    file_folder_func=lambda x: x\n",
    "    single_file=True\n",
    "f = PathSelector('.')\n",
    "display(f.accord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#Check if folder or file\n",
    "tiff_path=file_folder_func(f.accord.get_title(0))\n",
    "if not single_file:\n",
    "    #merge all tifs\n",
    "    cwd = os.getcwd()\n",
    "    tiffs_path = f'{tiff_path}/*_fcWM*.tif'\n",
    "    #tiffs=glob.glob(f\"{tiff_path}/{pattern}\")\n",
    "    #print(\"Processing the following files:\")\n",
    "    pprint.pprint(tiffs_path)\n",
    "    #combined_vrt=os.path.join(tiff_path, os.path.basename(tiff_path)+'.vrt')\n",
    "    get_ipython().system('gdalbuildvrt $tiff_path/water_mask_merged.vrt $tiffs_path')\n",
    "    combined_tif=os.path.join(tiff_path, os.path.basename(tiff_path)+'_merged_water_mask.tif')\n",
    "    #build_vrt(combined_vrt, tiffs)\n",
    "    #translate vrt to tif. There is a problem warping with the vrt.\n",
    "    cmd_translate=f\"gdal_translate -of GTiff {tiff_path}/water_mask_merged.vrt {combined_tif}\"\n",
    "    print(cmd_translate)\n",
    "    exitcode=os.WEXITSTATUS(os.system(cmd_translate))\n",
    "    if exitcode != 0:\n",
    "        print(\"Error in creating a mosaic from selected files\\n\")\n",
    "        print(f\"\\nCommand failed:\\n {cmd_translate}\")\n",
    "        assert exitcode == 0\n",
    "    tiff_path=combined_tif\n",
    "else:\n",
    "    pass # do nothing. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<br>\n",
    "<font face=\"Calibri\" size=\"4\"><b>Reproject tiffs from UTM to EPSG 4326:</b> </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(\"Choose your GDAL compatible precalculated HAND file using the file browser below:\")\n",
    "f = PathSelector('.')\n",
    "display(f.accord)        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#checking current coordinate reference system\n",
    "tiff_dir = os.path.dirname(tiff_path)\n",
    "info = (gdal.Info(tiff_path, options = ['-json']))\n",
    "info=info['coordinateSystem']['wkt']\n",
    "epsg = info.split('ID')[-1].split(',')[1].replace(']','')\n",
    "print(f\"EPSG code for Water Extent: {epsg}\")\n",
    "\n",
    "hand_dem=f.accord.get_title(0)\n",
    "info_hand = (gdal.Info(hand_dem, options = ['-json']))\n",
    "info_hand=info_hand['coordinateSystem']['wkt']\n",
    "epsg_hand = info_hand.split('ID')[-1].split(',')[1].replace(']','')\n",
    "print(f'EPSG for HAND: {epsg_hand}')\n",
    "\n",
    "\n",
    "#Reprojecting coordinate system\n",
    "if epsg != epsg_hand:\n",
    "    filename=tiff_path.split(tiff_dir)[1].split(\"/\")[1]\n",
    "    filenoext=os.path.splitext(filename)[0] #given vrt we want to force geotif output with tif extension\n",
    "    cmd_reproj=f\"gdalwarp -overwrite -t_srs EPSG:{epsg_hand} -r cubicspline -of GTiff {tiff_dir}/{filename} {tiff_dir}/reproj_{filenoext}.tif\"\n",
    "    print(cmd_reproj)\n",
    "    os.system(cmd_reproj)\n",
    "else:\n",
    "    filename=tiff_path.split(tiff_dir)[1].split(\"/\")[1]\n",
    "    filenoext=os.path.splitext(filename)[0] #given vrt we want to force geotif output with tif extension\n",
    "    if os.path.exists(f'{tiff_dir}/reproj_{filenoext}.tif'):\n",
    "        os.remove(f'{tiff_dir}/reproj_{filenoext}.tif')\n",
    "    os.symlink(f'{tiff_dir}/{filename}', f'{tiff_dir}/reproj_{filenoext}.tif')\n",
    "    \n",
    "# Building the virtual raster for Change Detection product(tiff)\n",
    "reprojected_flood_mask = f\"{tiff_dir}/reproj_{filenoext}.tif\"\n",
    "print(f\"Reprojected Flood Mask File: {reprojected_flood_mask}\")\n",
    "os.chdir(tiff_dir)\n",
    "\n",
    "pixels, lines=get_size(reprojected_flood_mask)\n",
    "print(f\"X-dimension: {pixels} Y-dimension: {lines}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#checking extent of the map\n",
    "info = (gdal.Info(reprojected_flood_mask, options = ['-json']))\n",
    "west,east,south,north=get_wesn(reprojected_flood_mask)\n",
    "print(f\"Retrieved Extent of Flood Extent (w/e/s/n):{west}, {east}, {south}, {north}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#Check if HAND is valid. \n",
    "hand_dem=f.accord.get_title(0)\n",
    "hand_dem_bb=bounding_box(hand_dem)\n",
    "if not bounding_box_inside_bounding_box(bounding_box(reprojected_flood_mask), hand_dem_bb):\n",
    "    print('Flood Extent Bounding Box:')\n",
    "    print(bounding_box(reprojected_flood_mask))\n",
    "    print('HAND bounding box:')\n",
    "    print(hand_dem_bb)\n",
    "    print('You can use BIG HAND Notebook to calculate HAND from a DEM.')\n",
    "    print('Image is not completely covered inside given HAND.')\n",
    "    print('If you continue your result may not be valid...')\n",
    "    if yesno(\"Do you want to continue?\"):\n",
    "        pass\n",
    "    else:\n",
    "        raise ValueError('Image is not completely covered inside given HAND.')\n",
    "\n",
    "#Clip HAND to the same size as the reprojected_flood_mask\n",
    "filename=os.path.basename(hand_dem)\n",
    "cmd_clip=f\"gdalwarp -overwrite -te {west} {south} {east} {north} -ts {pixels} {lines} -r lanczos  -of GTiff {hand_dem} {tiff_dir}/clip_{filename}\"\n",
    "print(cmd_clip)\n",
    "os.system(cmd_clip)\n",
    "hand_array=gis.readData(f\"{tiff_dir}/clip_{filename}\")\n",
    "if np.all(hand_array==0):\n",
    "    print('HAND is all zeros. HAND DEM does not cover the imaged area.')\n",
    "    raise ValueError # THIS SHOULD NEVER HAPPEN now that we are checking the bounding box. Unless HAND is bad.     \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "# Generating Flood Mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Pull Known Perennial Water Information from Public Repository"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\">All perennial Global Surface Water data is produced under the Copernicus Programme: Jean-Francois Pekel, Andrew Cottam, Noel Gorelick, Alan S. Belward, High-resolution mapping of global surface water and its long-term changes. Nature 540, 418-422 (2016). (doi:10.1038/nature20584). <b>We pull this layer to make sure all perennial water is accounted for in the surface water information that is used for Flood Depth Map calculation</b>.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#Get known Water Mask\n",
    "ths = known_water_threshold #30 #higher than 30% possibility \n",
    "known_water_mask=get_waterbody(reprojected_flood_mask, ths)\n",
    "if show_plots:\n",
    "    pl.matshow(known_water_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Grabbing Surface Water Extent Map Created in Lab 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\">Now we grab the Surface Water Extent Map that we created in Lab 2. </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#load and display change detection product from Hyp3\n",
    "hyp_map = gdal.Open(reprojected_flood_mask)\n",
    "change_map = hyp_map.ReadAsArray()\n",
    "\n",
    "#Initial mask layer generation\n",
    "for c in water_classes: # This allows more than a single water_class to be included in flood mask\n",
    "    change_map[change_map==c]=1\n",
    "\n",
    "mask=change_map==1\n",
    "flood_mask=np.bitwise_or(mask,known_water_mask) #add known water mask... #Added 20200921\n",
    "\n",
    "if show_plots:\n",
    "    pl.matshow(flood_mask);pl.title('Final Flood Mask')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Flood Depth Map Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<font face=\"Calibri\" size=\"3\"> Now we <b>add known water information to the SAR-derived surface water detection maps</b> and then we <b>generate our desired Flood Depth Product:</b> </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Calculate Flood Depth - Show Progress Bar\n",
    "flood_mask_labels, num_labels=ndimage.label(flood_mask)\n",
    "print(f'Detected {num_labels} water bodies...')\n",
    "object_slices=ndimage.find_objects(flood_mask_labels)\n",
    "if show_plots:\n",
    "    pl.matshow(flood_mask_labels);pl.colorbar();\n",
    "\n",
    "flood_depth=np.zeros(flood_mask.shape)\n",
    "print(f'Using estimator: {estimator}')\n",
    "for l in tqdm(range(1,num_labels)):#Skip first, largest label.\n",
    "    slices=object_slices[l-1] #osl label=1 is in object_slices[0]\n",
    "    min0=slices[0].start\n",
    "    max0=slices[0].stop\n",
    "    min1=slices[1].start\n",
    "    max1=slices[1].stop\n",
    "    flood_mask_labels_clip=flood_mask_labels[min0:max0, min1:max1] \n",
    "\n",
    "    flood_mask_clip=flood_mask[min0:max0, min1:max1].copy()\n",
    "    flood_mask_clip[flood_mask_labels_clip!=l]=0 #Maskout other flooded areas (labels)\n",
    "    hand_clip=hand_array[min0:max0, min1:max1] \n",
    "    if estimator.lower() == \"numpy\":\n",
    "        m=np.nanmean(hand_clip[flood_mask_labels_clip==l])\n",
    "        s=np.nanstd( hand_clip[flood_mask_labels_clip==l])\n",
    "        water_height=m+water_level_sigma*s\n",
    "    elif estimator.lower() == \"nmad\":\n",
    "        m=np.nanmean(hand_clip[flood_mask_labels_clip==l])\n",
    "        s=stats.median_abs_deviation(hand_clip[flood_mask_labels_clip==l], scale='normal', nan_policy='omit')\n",
    "        water_height=m+water_level_sigma*s\n",
    "    elif estimator.lower() == \"logstat\":\n",
    "        m=logstat(hand_clip[flood_mask_labels_clip==l], func=np.nanmean)\n",
    "        s=logstat(hand_clip[flood_mask_labels_clip==l])  \n",
    "        water_height=m+water_level_sigma*s\n",
    "    elif estimator.lower() == \"iterative\":\n",
    "        water_height=iterative(hand_clip, flood_mask_labels_clip==l, water_levels=iterative_bounds)            \n",
    "    else:\n",
    "        print(\"Unknown estimator selected for water height calculation.\")\n",
    "        raise ValueError\n",
    "    #if np.isnan(m) or np.isnan(s):\n",
    "    #    set_trace()\n",
    "    flood_depth_clip=flood_depth[min0:max0, min1:max1]\n",
    "    flood_depth_clip[flood_mask_labels_clip==l]=water_height-hand_clip[flood_mask_labels_clip==l]    \n",
    "\n",
    "#remove negative depths:\n",
    "flood_depth[flood_depth<0]=0\n",
    "if show_plots:\n",
    "    m=np.nanmean(flood_depth)\n",
    "    s=np.nanstd(flood_depth)\n",
    "    clim_min=max([m-2*s,0])\n",
    "    clim_max=min([m+2*s,5])\n",
    "    pl.matshow(flood_depth);pl.colorbar();pl.clim([clim_min,clim_max]);pl.title('Estimated Flood Depth')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "# Export Your Flood Depth Map as GeoTIFF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#Saving Estimated FD to geotiff\n",
    "gT = gis.getGeoTransform(f\"{tiff_dir}/clip_{filename}\")\n",
    "\n",
    "outfilename = tiff_path.split(tiff_dir)[1].split(\"/\")[1]\n",
    "srs_proj4=gdal_get_projection(f\"{tiff_dir}/clip_{filename}\")\n",
    "gis.writeTiff(flood_depth, gT, filename=\"_\".join(filter(None,[output_prefix,\"HAND_WaterDepth\",estimator,version,outfilename])), srs_proj4=srs_proj4, nodata=0, options = [\"TILED=YES\",\"COMPRESS=LZW\",\"INTERLEAVE=BAND\",\"BIGTIFF=YES\"])\n",
    "gis.writeTiff(flood_mask, gT, filename=\"_\".join(filter(None, [output_prefix,\"Flood_mask\",estimator,version,outfilename])), srs_proj4=srs_proj4, options = [\"TILED=YES\",\"COMPRESS=LZW\",\"INTERLEAVE=BAND\",\"BIGTIFF=YES\"])\n",
    "\n",
    "flood_mask[known_water_mask]=0\n",
    "flood_depth[np.bitwise_not(flood_mask)]=0\n",
    "gis.writeTiff(flood_depth, gT, filename=\"_\".join(filter(None,[output_prefix,\"HAND_FloodDepth\",estimator,version,outfilename])), nodata=0, srs_proj4=srs_proj4, options = [\"TILED=YES\",\"COMPRESS=LZW\",\"INTERLEAVE=BAND\",\"BIGTIFF=YES\"])\n",
    "print('Export complete.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Clean Up Temporary and Intermediate Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "outputs": [],
   "source": [
    "#clear some intermediate files\n",
    "try:\n",
    "    os.remove(reprojected_flood_mask)\n",
    "except:\n",
    "    pass\n",
    "try:\n",
    "    os.remove(f'{tiff_dir}/clip_{filename}')\n",
    "except:\n",
    "    pass\n",
    "try:\n",
    "    os.remove(f'{tiff_dir}/reproj_{filenoext}.tif')\n",
    "except:\n",
    "    pass\n",
    "try:\n",
    "    os.remove(f'{tiff_dir}/surface_water_map_clip.tif')\n",
    "except:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Version Log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "hideCode": false,
    "hidePrompt": false
   },
   "source": [
    "<font face=\"Calibri\" size=\"2\" color=\"gray\"> <i> Version 0.1.8 - Batu Osmanoglu, MinJeong Jo\n",
    "    \n",
    "<b>Change Log</b> <br>\n",
    "2021/01/24:<br>\n",
    "- Added `iterative` estimator. This method is based on `scipy.optimize.basinhopping` with bounds, which can be specified with the `iterative_bounds` parameter. It takes considerably longer to use iterative method as it tried to match the observed flood-extent pattern at different water levels.<br> \n",
    "2021/01/19:<br>\n",
    "- Minor cleanup and threshold implementation to `get_waterbody`. Also changed the dataset to 2019 (`downloads2019v2`). <br> \n",
    "2020/12/01:<br>\n",
    "- Added new statistical estimators for maximum water height: numpy, nmad or logstat. Numpy uses standard mean, and std.dev. functions. NMAD uses normalized mean absolute difference for sigma. See `scipy.stats.median_abs_deviation()` for details. logstat uses standard mean and std.dev functions after taking the logarithm of the data. See `logstat()` for details. <br>\n",
    "2020/11/09: <br>\n",
    "-Changed known water source, Occurrence 2019 vertion. Added a threshold variable. <br>\n",
    "2020/10/22: <br>\n",
    "-BugFix: During reproj EPSG code was wrongly identified. Also if EPSG code is read-wrong the -s_srs flag in gdal_warp causing the reprojection to fail. Fixed both. <br>\n",
    "-Testing: Replaced standard deviation with normalized mean absolute deviation. For large riverine floods, large water objects result in overestimation of sigma resulting in deeper than expected water depths.  <br>\n",
    "-Feat: Removing negative water depths in the final map. <br>\n",
    "2020/10/10: <br>\n",
    "-BugFix: Looks like with the recent updates `gdal.Info(tiff_path, options = ['-json'])` returns a dict instead of a string. Fixed the collection of projection based on this. <br>\n",
    "-Feat: Allowing to continue even if HAND is smaller than image. This is useful if SAR image covers significant amount of ocean etc. <br>\n",
    "-BugFix: gis.readData() was failing to read the VRT generated in get_waterbody. gdalwarp outputs a GeoTif now. <br>\n",
    "2020/10/01: <br>\n",
    "- Feat: Moving away from repetitive argwhere to ndimage.find_objects results in ~6000x faster calculation.\n",
    "2020/09/30: <br>\n",
    "- BugFix: The known water mask and input mask was merged using a sum operator instead of a bitwise_or operator.<br> \n",
    "2020/09/20:<br>\n",
    "- BugFix: Added known water download and addition to the mask. This helps to make sure known water bodies are handled as a single object. Removed morphological filters also. <br>\n",
    "2020/09/17:<br>\n",
    "- First version. \n",
    "</i></font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_code_all_hidden": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.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": 2
}