Big_Hand_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false
   },
   "source": [
    "<img src=\"https://radar.community.uaf.edu/wp-content/uploads/sites/667/2021/03/HydroSARbanner.jpg\" width=\"100%\" />\n",
    "<font size=\"7\"> <b> Calculating HAND from DEM</b></font>\n",
    "\n",
    "<font size=\"5\">  How to obtain Height Above Nearest Drainage using a BIG Digital Elevation Model with PySheds </font>\n",
    "\n",
    "<br>\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.13 - 2022/02/25 <br>\n",
    "<b>Change Log</b><br>\n",
    "See bottom of the notebook.<br>\n",
    "</font> \n",
    "<font color='rgba(0,0,200,0.2)'> <b>Contact: </b> batuhan.osmanoglu@nasa.gov </font>\n",
    "\n",
    "\n",
    "<hr>\n",
    "<font face=\"Calibri\">\n",
    "\n",
    "<font size=\"5\"> <b> 0. Importing Relevant Python Packages </b> </font>\n",
    "\n",
    "<font size=\"3\"> The first step in any notebook is to import the required Python libraries into the Jupyter environment. In this notebooks we use the following libraries:\n",
    "<ol type=\"1\">\n",
    "    <li> <b><a href=\"https://www.gdal.org/\" target=\"_blank\">GDAL</a></b> is a software library for reading and writing raster and vector geospatial data formats. It includes a collection of programs tailored for geospatial data processing. Most modern GIS systems (such as ArcGIS or QGIS) use GDAL in the background.</li>\n",
    "    <li> <b><a href=\"http://www.numpy.org/\" target=\"_blank\">NumPy</a></b> is one of the principal packages for scientific applications of Python. It is intended for processing large multidimensional arrays and matrices, and an extensive collection of high-level mathematical functions and implemented methods makes it possible to perform various operations with these objects. </li>\n",
    "    <li> <b><a href=\"https://www.tutorialspoint.com/matplotlib/matplotlib_pylab_module.htm\" target=\"_blank\">PyLab</a></b> is a procedural interface to the Matplotlib object-oriented plotting library. Matplotlib is the whole package; matplotlib.pyplot is a module in Matplotlib; and PyLab is a module that gets installed alongside Matplotlib.</li>\n",
    "    <li> <b><a href=\"https://github.com/mdbartos/pysheds\" target=\"_blank\">PySHEDS</a></b> is a new python module developed by Matt Bartos with the motto \"Simple and fast watershed delineation in python.\" We use the Height Above Nearest Drainage (HAND) implementation in PySHEDS.</li>\n",
    "    <li> <b><a href=\"https://github.com/Toblerity/Fiona\" target=\"_blank\"> Fiona </a></b> reads and writes geographic data files and thereby helps Python programmers integrate geographic information systems with other computer systems. Fiona contains extension modules that link the Geospatial Data Abstraction Library (GDAL). </li>\n",
    "    <li> <b><a href=\"https://github.com/Toblerity/Shapely\" target=\"_blank\"> Shapely </a></b>allows for manipulation and analysis of geometric objects in the Cartesian plane. Shapely is not concerned with data formats or coordinate systems, but can be readily integrated with packages that are.</li>\n",
    "    <li> <b><a href=\"https://geopandas.org/\" target=\"_blank\"> GeoPandas </a></b> is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types. Geometric operations are performed by shapely. Geopandas further depends on fiona for file access and descartes and matplotlib for plotting.</li>\n",
    "    <li> <b><a href=\"https://www.astropy.org/\" target=\"_blank\"> The Astropy Project </a></b> is a community effort to develop a core package for astronomy using the Python programming language and improve usability, interoperability, and collaboration between astronomy Python packages. The core astropy package contains functionality aimed at professional astronomers and astrophysicists, but may be useful to anyone developing astronomy software.</li>\n",
    "    <li> <b><a href=\"https://github.com/tqdm/tqdm\" target=\"_blank\"> tqdm </a></b> is a smart progress meter that allows easy addition of a loop counter.</li>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Setup Environment\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "from osgeo import gdal,osr\n",
    "import pylab as pl\n",
    "from functools import partial\n",
    "import warnings #Suppress warnings on occasion\n",
    "import tempfile #for creation of temporary folder\n",
    "import urllib   #for data download\n",
    "import zipfile  #zipfile\n",
    "import time     #to time HAND nan-fill duration\n",
    "\n",
    "from IPython.core.debugger import set_trace #Enable if you like to debug and add set_trace() where you want debugger\n",
    "\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",
    "#This line below is for direct download\n",
    "from IPython.display import FileLink\n",
    "#This is to enable interactive plotting in Jupyterlab\n",
    "try:\n",
    "    import ipympl\n",
    "except:\n",
    "    !pip install ipympl\n",
    "    import ipympl\n",
    "#All packages below are used to support pysheds. \n",
    "try:\n",
    "    import pysheds\n",
    "except:\n",
    "    !pip install pysheds\n",
    "    import pysheds\n",
    "from pysheds.grid import Grid\n",
    "from pysheds.pgrid import Grid as Pgrid\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",
    "    with warnings.catch_warnings():\n",
    "        warnings.filterwarnings(\"ignore\")#, category=FutureWarning)    \n",
    "        import pyproj\n",
    "except:\n",
    "    !pip install pyproj\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.filterwarnings(\"ignore\")#, category=FutureWarning)\n",
    "        import pyproj #Pyproj is showing a lot of future warnings. \n",
    "        #/opt/conda/lib/python3.7/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n",
    "#This is for masking the DEM based on hydrobasins shapefile.\n",
    "try:\n",
    "    import fiona\n",
    "except:\n",
    "    !pip install fiona\n",
    "    import fiona\n",
    "import fiona.crs\n",
    "try:\n",
    "    import shapely\n",
    "except:\n",
    "    !pip install shapely\n",
    "    import shapely\n",
    "from shapely.ops import transform\n",
    "try:\n",
    "    import geopandas as gpd\n",
    "except:\n",
    "    !pip install geopandas\n",
    "    import geopandas as gpd\n",
    "#Astropy and scipy ndimage is to interpolate over nans in final HAND. \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",
    "#pretty loop counter\n",
    "try:\n",
    "    from tqdm.auto import tqdm\n",
    "except:\n",
    "    !pip install tqdm\n",
    "    from tqdm.auto import tqdm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=\"5\"> <b> 1. Define convenience functions </b> </font>\n",
    "\n",
    "<font size=\"3\"> Here we define some functions for later convenience.\n",
    "    \n",
    "<ol type=\"1\">\n",
    "    <li> <b>fiona_read_vectorfile</b> returns a list of shapes (and optionally properties) using fiona. </li>\n",
    "    <li> <b>fiona_write_vectorfile</b> write a vectorfile containing shapes (and optionally properties) with fiona. </li>\n",
    "    <li> <b>intersect</b> returns polygons from multiple 'geometries' read by fiona.</li>\n",
    "    <li> <b>reproject</b> reprojects a given vector file to another coordinate reference system (CRS). </li>\n",
    "    <li> <b>transform_shape</b> transforms a single geometry to another coordinate reference system (CRS).</li>\n",
    "    <li> <b>xy2coord</b> converts pixel index to position based on geotransform. </li>     \n",
    "    <li> <b>get_projection</b> returns the spatial reference system in wkt, proj4, or epsg formats.</li>\n",
    "    <li> <b>gdal_get_geotransform</b> returns the geotransform of the dataset using GDAL. </li>\n",
    "    <li> <b>gdal_get_size</b> returns width and height for a given raster file. </li>\n",
    "    <li> <b>gdal_bounding_box</b> returns the bounding box in geocoded coordinates for a given raster file. </li>    \n",
    "    <li> <b>gdal_write</b> is used to export the generated HAND product using GDAL. </li>\n",
    "    <li> <b>fill_nan_based_on_DEM</b> is used to fill small not-a-number areas in final HAND product if necessary. </li>\n",
    "    <li> <b>calculate_hand</b> returns the height above nearest drainage for a given DEM using the pySheds library. </li>\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     1,
     20,
     30,
     63,
     78,
     87,
     94,
     103,
     127,
     135,
     144,
     152,
     188
    ]
   },
   "outputs": [],
   "source": [
    "# Define convenience functions\n",
    "def fiona_read_vectorfile(vectorfile, get_property=None):\n",
    "    \"\"\"shapes=fiona_read_vectorfile(vectorfile, get_property=None)\n",
    "       shapes, props=fiona_read_vectorfile(vectorfile, get_property='Property_Name')\n",
    "       Returns a list of shapes (and optionally properties) using fiona.\n",
    "       \n",
    "       vectorfile: any fiona compatible vector file. \n",
    "       get_property: String for the property to be read. \n",
    "       shapes: List of vector \"geometry\"\n",
    "       props:  List of vector \"properties\"\n",
    "    \"\"\"\n",
    "    with fiona.open(vectorfile, \"r\") as shpf:\n",
    "        shapes   = [ feature[\"geometry\"] for feature in shpf ]\n",
    "        print(f\"Number of shapes loaded: {len(shapes)}\")\n",
    "        if get_property is not None:\n",
    "            props = [ feature[\"properties\"][get_property] for feature in shpf ]\n",
    "            return shapes, props\n",
    "        else:\n",
    "            return shapes    \n",
    "\n",
    "def fiona_write_vectorfile(shapes, vectorfile, crs=fiona.crs.from_epsg(4326), driver='ESRI Shapefile', schema_type='Polygon'):\n",
    "    if schema_type=='Polygon':\n",
    "        schema={'geometry': 'Polygon',\n",
    "                'properties': {}}        \n",
    "    with fiona.open(vectorfile, 'w',crs=crs,driver=driver, schema=schema) as output:\n",
    "        for s in shapes:\n",
    "            if schema_type=='Polygon':\n",
    "                sp= shapely.geometry.Polygon(s)\n",
    "            output.write({'geometry':shapely.geometry.mapping(sp),'properties': {}})    \n",
    "        \n",
    "def intersect(shapes, polygon, properties=None):\n",
    "    \"\"\"\n",
    "    polygons=intersect(shapes, polygon, properties=None)\n",
    "    Returns polygons from multiple 'geometries' read by fiona.\n",
    "    \n",
    "    shapes: shapes returned by fiona_read_vectorfile()\n",
    "    polygon: a single polygon to intersect with shapes\n",
    "    properties: If not none, returns the property value instead of polygon geometry.\n",
    "    \"\"\"\n",
    "    #first loop to split multi polygons to single polygons\n",
    "    polygons=[]\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter(\"ignore\")    \n",
    "        for k,shape in enumerate(tqdm(shapes)):\n",
    "            if shape['type']=='MultiPolygon':\n",
    "                for l,p in enumerate(shape['coordinates']):\n",
    "                    s=shapely.geometry.Polygon(p[0])\n",
    "                    if polygon.intersects(s) and properties is None:            \n",
    "                        polygons.append(s)\n",
    "                    elif polygon.intersects(s) and properties is not None:\n",
    "                        if np.isscalar(properties[k]):\n",
    "                            polygons.append(properties[k])\n",
    "                        else:\n",
    "                            polygons.append(properties[k][l])\n",
    "                    \n",
    "            elif shape['type']=='Polygon':\n",
    "                s=shapely.geometry.Polygon(shape['coordinates'][0])\n",
    "                if polygon.intersects(s) and properties is None:\n",
    "                    polygons.append(s)\n",
    "                elif polygon.intersects(s) and properties is not None:\n",
    "                    polygons.append(properties[k])\n",
    "    return polygons\n",
    "\n",
    "def reproject(vector_file, output_crs, output_file=None):\n",
    "    \"\"\"\n",
    "    output_file=reproject(vector_file, output_crs, output_file=None)\n",
    "    Reprojects a given vector file to another reference frame (CRS). \n",
    "    vector_file: Any vector file that can be opened with GeoPandas\n",
    "    output_crs: A rasterio opened crs (e.g. dem.crs)\n",
    "    output_file: if not defined, defaults to vector_file[:-4]+'_warp.shp'. \n",
    "    \"\"\"\n",
    "    v=gpd.GeoDataFrame.from_file(vector_file)\n",
    "    warp=v.to_crs(output_crs)\n",
    "    if output_file is None:\n",
    "        output_file=vector_file[:-4]+'_warp.shp'\n",
    "    warp.to_file(output_file)\n",
    "    return output_file\n",
    "\n",
    "def transform_polygon(polygon, s_srs='EPSG:4269', t_srs='EPSG:4326'):\n",
    "    shp_geom = shapely.geometry.Polygon(polygon)    \n",
    "    project = pyproj.Transformer.from_proj(\n",
    "        pyproj.Proj(init=s_srs), # source coordinate system\n",
    "        pyproj.Proj(init=t_srs)) # destination coordinate system\n",
    "    \n",
    "    # polygon is a shapley Polygon\n",
    "    return transform(project.transform, shp_geom)  # apply projection\n",
    "\n",
    "def transform_shape(shape, s_srs='epsg:4326', t_srs='epsg:4326'):\n",
    "    transformation=partial(\n",
    "               pyproj.transform,\n",
    "               pyproj.Proj(init=s_srs), #source coordinate system\n",
    "               pyproj.Proj(init=t_srs)) #destination coordinate system\n",
    "    return shapely.ops.transform(transformation, shape)\n",
    "\n",
    "def xy2coord(x,y,gT):\n",
    "    '''\n",
    "    lon,lat=xy2coord(x,y,geoTransform)\n",
    "    converts 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 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",
    "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 gdal_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",
    "\n",
    "def gdal_bounding_box(filename):\n",
    "    \"\"\"\n",
    "    ((lon1,lat1), (lon2,lat2), (lon3,lat3), (lon4,lat4))=bounding_box('/path/to/file')\n",
    "    \"\"\"\n",
    "    gT=gdal_get_geotransform(filename)\n",
    "    width, height=gdal_get_size(filename)     \n",
    "    return (xy2coord(0,0,gT), xy2coord(width,0,gT), xy2coord(width, height,gT), xy2coord(0, height,gT))\n",
    "\n",
    "def gdal_write(ary, geoTransform, fileformat=\"GTiff\", filename='jupyter_rocks.tif', data_format=gdal.GDT_Float64, nodata=None, srs_proj4='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs', options = [\"TILED=YES\",\"COMPRESS=LZW\",\"INTERLEAVE=BAND\",\"BIGTIFF=YES\"], build_overviews=True):\n",
    "    '''gdal_write(ary, geoTransform, format=\"GTiff\", filename='jupyter_rocks.tif', data_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, data_format, options)\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",
    "    if build_overviews:\n",
    "        ds.BuildOverviews(\"NEAREST\", [2, 4, 8, 16, 32, 64, 128, 256])\n",
    "    ds = None\n",
    "    print(\"File written to: \" + filename);\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 fill_nan_based_on_DEM(arr, dem):\n",
    "    \"\"\"\n",
    "    filled_arr=fill_nan_based_on_DEM(arr, dem)\n",
    "    Fills Not-a-number values in arr using astropy. \n",
    "    \"\"\"    \n",
    "    hond = dem - arr; #height of nearest drainage \n",
    "    kernel = astropy.convolution.Gaussian2DKernel(x_stddev=3) #kernel x_size=8*stddev\n",
    "    arr_type=hond.dtype          \n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter(\"ignore\")\n",
    "        while np.any(np.isnan(hond)):\n",
    "            hond = astropy.convolution.interpolate_replace_nans(hond.astype(float), kernel, convolve=astropy.convolution.convolve)\n",
    "    my_mask=np.isnan(arr)\n",
    "    arr[my_mask]=dem[my_mask]-hond[my_mask]\n",
    "    return arr.astype(arr_type) \n",
    "\n",
    "def calculate_hand(dem, dem_gT, dem_proj4, mask=None, verbose=False, acc_thresh=100):\n",
    "    \"\"\"\n",
    "    hand=calculate_hand(dem, dem_gT, dem_proj4, mask=None, verbose=False)\n",
    "    Calculate the height above nearest drainage using pySHEDS library. This is done over a few steps:\n",
    "\n",
    "    Fill_Depressions fills depressions in a DEM (regions of cells lower than their surrounding neighbors).\n",
    "    Resolve_Flats resolves drainable flats in a DEM.\n",
    "    FlowDir converts the DEM to flow direction based on dirmap.\n",
    "    Accumulation converts from flow direction to flow accumulation.\n",
    "    Compute_Hand is used to convert directions to height above nearest drainage.\n",
    "    \n",
    "    NaN values are filled at the end of resolve_flats and final steps. \n",
    "    \n",
    "    Inputs:\n",
    "      dem=Numpy array of Digital Elevation Model (DEM) to convert to HAND. \n",
    "      dem_gT= GeoTransform of the input DEM\n",
    "      dem_proj4=Proj4 string of DEM\n",
    "      mask=If provided parts of DEM can be masked out. If not entire DEM is evaluated. \n",
    "      verbose=If True, provides information about where NaN values are encountered. \n",
    "      acc_thresh=Accumulation threshold. By default is set to 100. If none, \n",
    "                 mean value of accumulation array (acc.mean()) is used. \n",
    "    \"\"\"\n",
    "\n",
    "    #Specify  directional mapping\n",
    "             #N , NE , E ,SE,S,SW, W , NW\n",
    "    dirmap = (64, 128, 1, 2, 4, 8, 16, 32) \n",
    "    #Load DEM into pySheds\n",
    "    if type(dem_gT)==Affine:\n",
    "        aff=dem_gT\n",
    "    else:\n",
    "        aff=Affine.from_gdal(*tuple(dem_gT))\n",
    "    if mask is None:\n",
    "        mask=np.ones(dem.shape, dtype=np.bool)\n",
    "    grid=Pgrid() #(shape=dem.shape,affine=aff, crs=dem_proj4, mask=mask)\n",
    "    grid.add_gridded_data(dem, data_name='dem',affine=aff, crs=dem_proj4, mask=mask)        \n",
    "    #Fill Depressions\n",
    "    grid.fill_depressions('dem', out_name='flooded_dem')\n",
    "    if verbose:\n",
    "        set_trace()\n",
    "    if np.any(np.isnan(grid.flooded_dem)):\n",
    "        if verbose:\n",
    "            print('NaN:fill_depressions')\n",
    "            grid.flooded_dem=fill_nan(grid.flooded_dem)     \n",
    "    #Resolve_Flats \n",
    "    #Please note that Resolve_Flats currently has an open bug and can fail on occasion. https://github.com/mdbartos/pysheds/issues/118\n",
    "    try:\n",
    "        grid.resolve_flats('flooded_dem', out_name='inflated_dem')\n",
    "    except:\n",
    "        grid.inflated_dem=grid.flooded_dem\n",
    "    #if np.sum(np.isnan(grid.inflated_dem))<dem.size*0.5: #if nans account for less than 50% of the dem nanfill. \n",
    "    #    if verbose:\n",
    "    #        print('NaN:resolve_flats but less than 50%. Applying large value')\n",
    "    #    grid.inflated_dem=fill_nan(grid.inflated_dem)\n",
    "    if np.any(np.isnan(grid.inflated_dem)):\n",
    "        if verbose:\n",
    "            print('NaN:resolve_flats replacing with inflated_dem')\n",
    "        #grid.inflated_dem=fill_nan(grid.inflated_dem)                    \n",
    "        grid.inflated_dem[np.isnan(grid.inflated_dem)] = dem[np.isnan(grid.inflated_dem)]#10000  # setting nan to 10.000 to ensure drainage\n",
    "        ### Ref: https://github.com/mdbartos/pysheds/issues/90\n",
    "    #Obtain flow direction\n",
    "    grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap, apply_mask=True)\n",
    "    if np.any(np.isnan(grid.dir)):\n",
    "        if verbose:\n",
    "            print('NaN:flowdir')\n",
    "            grid.dir=fill_nan(grid.dir)     \n",
    "    #Obtain accumulation\n",
    "    grid.accumulation(data='dir', dirmap=dirmap, out_name='acc')\n",
    "    if np.any(np.isnan(grid.acc)):\n",
    "        if verbose:\n",
    "            print('NaN:accumulation')\n",
    "            grid.acc=fill_nan(grid.acc)     \n",
    "    #Generate HAND\n",
    "    if acc_thresh is None:\n",
    "        acc_thresh=grid.acc.mean()\n",
    "    #grid.compute_hand('dir', 'inflated_dem', grid.acc >100, out_name='hand')\n",
    "    #Copy HAND as an array. \n",
    "    #hand=grid.view('hand')    \n",
    "    hand = grid.compute_hand('dir', 'inflated_dem', grid.acc > acc_thresh, inplace=False)\n",
    "    if np.any(np.isnan(hand)):\n",
    "        if verbose:\n",
    "            print('NaN:compute_hand')\n",
    "            #attempt to fill low-lying flat areas with zeros. In radar DEMs vegetation alongside river, can trap\n",
    "            #the river and not let any water go into the river. This was seen in Bangladesh with SRTM 1 arcsec\n",
    "            #and NASADEM at Hydro Basin with ID: 4120928640\n",
    "            \n",
    "            #get nans inside masked area and find mean height for pixels outside the nans (but inside basin mask)\n",
    "            valid_nanmask=np.logical_and(mask, np.isnan(hand))\n",
    "            valid_mask   =np.logical_and(mask, ~np.isnan(hand)) \n",
    "            mean_height=grid.inflated_dem[valid_mask].mean()\n",
    "            #calculate gradient and set mean gradient magnitude as threshold for flatness. \n",
    "            g0,g1=np.gradient(grid.inflated_dem);\n",
    "            gMag=np.sqrt(g0**2+g1**2)\n",
    "            gMagTh=np.min(1, np.mean(gMag*np.isnan(hand)) ) # Make sure this threshold is not too high. We don't want to set rough surfaces to zero.\n",
    "            \n",
    "            #define low lying (<mean) pixels inside valid area. \n",
    "            #valid_flats=np.logical_and(valid_nanmask, grid.dir==0) #I thought grid.dir=0 meant flats. But this is not the case always apparently. \n",
    "            valid_flats=np.logical_and(valid_nanmask, gMag<gMagTh)\n",
    "            valid_low_flats=np.logical_and(valid_flats, grid.inflated_dem<mean_height)            \n",
    "            hand[ valid_low_flats ]=0\n",
    "        if np.any(np.isnan(hand)):\n",
    "            grid.hand=fill_nan(hand) \n",
    "    return hand\n",
    "\n",
    "def point_coordinates_to_geometry(coordinates, geometry_type='Polygon'):\n",
    "    if geometry_type.lower() == 'polygon':\n",
    "      return shapely.geometry.Polygon(coordinates)\n",
    "    else:\n",
    "      raise NotImplementedError\n",
    "def print_duration(start):\n",
    "    print('   - Processing time: %6.2f minutes' % ((time.time() - start)/60.0))\n",
    "    print('-----------------------------------------------------------------------------------')\n",
    "    print(' ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=\"5\"> <b> 1.1 Path Selector </b> </font>\n",
    "\n",
    "<font size=\"3\"> Path selector is used to select a raster DEM as as input to the HAND calculation. \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "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": {},
   "source": [
    "<font size=\"5\"> <b> 2. Load DEM & Set processing parameters </b> </font>\n",
    "\n",
    "<font size=\"3\">Load and display the digital elevation model. The DEM can be uploaded using the file Jupyterhub File Browser before selection below. \n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define some processing parameters\n",
    "version=\"0.1.13\" #used in filenames as a suffix\n",
    "show_plots=True #set True if you like to see plots providing more information on steps.\n",
    "temporary_folder_object=tempfile.TemporaryDirectory()\n",
    "temp_dir=temporary_folder_object.name #Folder name to be used in generating temporary files\n",
    "hybas_dir='~/external_data' # if you do not want to keep any hybas files set it as hybas_dir=temp_dir\n",
    "gshhg_dir='~/external_data' #if you do not want to keep the coastline file, set it as gshhg_dir=temp_dir\n",
    "debug=False #If true print more detailed messages. Turn false when using IPDB to debug. Otherwise plots won't be visible within IPDB. \n",
    "accumulation_threshold=None # This sets how large of an accumulation area is used for HAND. If too small, we get a very fine river network, which can be noisy. If too high, we get a very smooth HAND... \n",
    "                           # Recommended values None (for automatic) or 100. \n",
    "pad_width=1 # Padding applied to the hydrobasins polygons for HAND processing. At least 1 pixel is recommended. \n",
    "nodata_fill_value=np.finfo(float).eps # If set to None, nodata/ocean values are not touched (likely returning np.nan). If a value is specified, it will be used to fill oceans and no-data areas in input dem. \n",
    "dem_nodata_value=0 # If set to None, will read from dem.nodatavals property. If set, uses that value as a mask to skip nan-interpolation \n",
    "# define URLs for external data used in this project. \n",
    "gshhg_url='http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-shp-2.3.7.zip'\n",
    "hybas_extents_url='https://www.dropbox.com/s/fthjqjnxj829d7p/hybas_extent_v1c.gpkg?dl=1'\n",
    "hybas_links={'af':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AABSBGFylsZ9KoG8zYRvOTzqa/HydroBASINS/standard/af/hybas_af_lev12_v1c.zip?dl=1',\n",
    "             'eu':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AADULrBSkGy5dHOZ8vMxWpWxa/HydroBASINS/standard/eu/hybas_eu_lev12_v1c.zip?dl=1',\n",
    "             'si':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AABtI2KbgItfLp4jmHcvZhDea/HydroBASINS/standard/si/hybas_si_lev12_v1c.zip?dl=1',\n",
    "             'as':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AADWZKiGaCncO5JdRLmkIduMa/HydroBASINS/standard/as/hybas_as_lev12_v1c.zip?dl=1',\n",
    "             'au':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AAA5lwuZZ5EZsxrx_EBQGW3ma/HydroBASINS/standard/au/hybas_au_lev12_v1c.zip?dl=1',\n",
    "             'sa':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AABPzWxd07pmshjZl6Y0NPXNa/HydroBASINS/standard/sa/hybas_sa_lev12_v1c.zip?dl=1',\n",
    "             'na':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AAA1ofV7PhSY_x7vQluubYyNa/HydroBASINS/standard/na/hybas_na_lev12_v1c.zip?dl=1',\n",
    "             'ar':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AADaA0icxaPYgaQGuLbSaKfna/HydroBASINS/standard/ar/hybas_ar_lev12_v1c.zip?dl=1',\n",
    "             'gr':'https://www.dropbox.com/sh/hmpwobbz9qixxpe/AACNOTXj-M1T-rpz5k_QJd6Ka/HydroBASINS/standard/gr/hybas_gr_lev12_v1c.zip?dl=1'}\n",
    "if show_plots:\n",
    "    %matplotlib inline \n",
    "    # set to inline to revert back to Jupyter default. Jupyterlab does not support interactive out of the box. \n",
    "    #https://stackoverflow.com/questions/51922480/javascript-error-ipython-is-not-defined-in-jupyterlab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<font face=\"Calibri\" size=\"4\"><b>2.1 Select your DEM using the file-tree below:</b> </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Display file selector to select DEM\n",
    "print(\"Choose your GDAL compatible DEM using the file browser below:\")\n",
    "f = PathSelector('.')\n",
    "display(f.accord)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b>2.2 Obtain DEM Parameters</b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\">We obtain several parameters for the DEM like it's geotransform, projection and bounding box. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#Obtain DEM parameters like, projection, geoTransform, bounding box etc. \n",
    "t=time.time(); t0=t\n",
    "dem_file=f.accord.get_title(0)\n",
    "\n",
    "try:\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
    "        #load and display dem\n",
    "        print(f\"Selected DEM: {dem_file}\")\n",
    "        dem_gT=gdal_get_geotransform(dem_file)\n",
    "        dem_proj4=get_projection(dem_file)\n",
    "\n",
    "        import rasterio.mask\n",
    "        dem=rasterio.open(dem_file)\n",
    "        bb=gdal_bounding_box(dem_file)\n",
    "        bb_poly=list(bb)\n",
    "        bb_poly.append(bb[0])\n",
    "        dem_poly=shapely.geometry.Polygon(bb_poly)    \n",
    "        dem_poly_wgs84=transform_shape(dem_poly, s_srs=dem.crs.to_string())\n",
    "except:\n",
    "    # don't go beyond here with Run All if above did not work. \n",
    "    print('Please select a GDAL compatible DEM using the file-tree above.')\n",
    "    assert False\n",
    "print_duration(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"5\"> <b>3. Load Basins</b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\">Here we load the apropriate basin information.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Download Basin Extents (hybas_extent_v1c.gpkg) from Batu's dropbox\n",
    "t=time.time()\n",
    "#Expand Tilda in folder names\n",
    "temp_dir=os.path.expanduser(temp_dir)\n",
    "hybas_dir=os.path.expanduser(hybas_dir)\n",
    "gshhg_dir=os.path.expanduser(gshhg_dir)\n",
    "if not os.path.exists(hybas_dir):\n",
    "    os.mkdir(hybas_dir)\n",
    "if not os.path.exists(gshhg_dir):\n",
    "    os.mkdir(gshhg_dir)    \n",
    "extent_file=os.path.join(hybas_dir,'hybas_extent_v1c.gpkg')\n",
    "if not os.path.exists(extent_file):\n",
    "    #!wget -O hybas_extent_v1c.gpkg https://www.dropbox.com/s/fthjqjnxj829d7p/hybas_extent_v1c.gpkg?dl=1\n",
    "    urllib.request.urlretrieve(hybas_extents_url, extent_file)\n",
    "    \n",
    "#find which hydrobasin data to download. \n",
    "pf_dict={1:'af',2:'eu',3:'si',4:'as',5:'au',6:'sa',7:'na',8:'ar',9:'gr'} #PF=Pfafstetter Code, https://www.hydrosheds.org/images/inpages/HydroBASINS_TechDoc_v1c.pdf\n",
    "pf_desc={'af':'Africa', 'eu':'Europe', 'si':'Siberia', 'as':'Asia', 'au':'Australia', 'sa':'South America', 'na':'North America', 'ar':'Arctic', 'gr':'Greenland'}\n",
    "\n",
    "#read extent shapes\n",
    "shapes,pf_codes=fiona_read_vectorfile(extent_file, get_property='PF_CODE')\n",
    "\n",
    "#find intersecting shapes\n",
    "polygons=intersect(shapes, dem_poly_wgs84, properties=pf_codes)\n",
    "#Find the correct Pfafstetter code \n",
    "if any(polygons):\n",
    "    if len(np.unique(polygons))==1: #polygons==[polygons[0]]:\n",
    "        pf_str=pf_dict[polygons[0]]\n",
    "        print(f'Detected hydrobasin location: {pf_desc[pf_str]}')\n",
    "    else:\n",
    "        print(f'The DEM is intersecting with {len(np.unique(polygons))} continents.')\n",
    "        print('Please select the continent you would like to process:')\n",
    "        print(pf_dict)        \n",
    "        pf_str=input()\n",
    "else: \n",
    "    print(\"Can not find a compatible hydrobasins area for this DEM. If it's a small DEM, try regular HAND instead of Big HAND.\")\n",
    "    assert False    \n",
    "#Download Hydrobasins data\n",
    "hybas_zipfile=os.path.join(hybas_dir,f\"hybas_{pf_str}_lev12_v1c.zip\")\n",
    "hybas_file=os.path.join(hybas_dir, f\"hybas_{pf_str}_lev12_v1c.shp\")\n",
    "if not os.path.exists(hybas_zipfile) and not os.path.exists(hybas_file):\n",
    "    #!wget -O {hybas_zipfile} {hybas_links[pf_str]}\n",
    "    urllib.request.urlretrieve(hybas_links[pf_str], hybas_zipfile)\n",
    "if not os.path.exists(hybas_file): \n",
    "    #!unzip -o {hybas_zipfile} # The HydroBASINS_TechDoc_v1c.pdf is common and requires overwriting.\n",
    "    with zipfile.ZipFile(hybas_zipfile, 'r') as zip_ref:          \n",
    "        zip_ref.extractall(path=hybas_dir) #overwrites HydroBASINS_TechDoc_v1c.pdf by default\n",
    "    os.remove(hybas_zipfile)\n",
    "print_duration(t)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b>3.1 Match basins to DEM projection</b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\"> Use `reproject()` function to convert the hydrobasin file (`hybas_file`) to the same projection as the DEM if necessary.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Convert the Hydrobasins to the DEM projection. \n",
    "t=time.time()\n",
    "if dem.crs.to_string() != 'EPSG:4326':\n",
    "    print(\"DEM and Hydrobasins projections differ.\")\n",
    "    print(\"Reprojecting Hydrobasins shapefile.\")\n",
    "    if not os.path.exists(hybas_file[:-4]+'_epsg4326.shp'):\n",
    "        hybas_file=reproject(hybas_file, dem.crs, output_file=hybas_file[:-4]+'_epsg4326.shp')\n",
    "    else:\n",
    "        hybas_file=hybas_file[:-4]+'_epsg4326.shp'  \n",
    "    print(f\"Output File: {hybas_file}\")\n",
    "print_duration(t)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b>3.2 Read and Intersect HydroBasins</b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\">Read all polygons in the HydroBasins file and intersect with the DEM bounding box stored in `dem_poly`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Find basins intersecting the DEM.\n",
    "t=time.time()\n",
    "shapes,hybas_id=fiona_read_vectorfile(hybas_file, get_property='HYBAS_ID')\n",
    "polygons=intersect(shapes, dem_poly)\n",
    "polygon_ids=intersect(shapes, dem_poly, properties=hybas_id)\n",
    "print(f\"Number of polygons intersecting the DEM: {len(polygons)}\")\n",
    "print_duration(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b>3.3 Calculate HAND</b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\">Loop over the intersecting `polygons` and calculate hand. We are using Python's assignment by reference to make place the hand values for a small area `h` in large `hand` array. \n",
    "    \n",
    "Ref: https://docs.python.org/2.0/ref/assignment.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Loop over each basin and calculate HAND\n",
    "t=time.time()\n",
    "if dem_nodata_value is None:\n",
    "    dem_nodata_value=dem.nodatavals[0]\n",
    "if dem_nodata_value is None:\n",
    "    print('DEM does not have a defined no-data value.')\n",
    "    print('Assuming all valid pixels. If not, expect long processing times.')\n",
    "else:\n",
    "    dem_nodata_mask=dem.read(1)==dem_nodata_value\n",
    "    \n",
    "hand=np.zeros(dem.shape)\n",
    "hand[:]=np.nan #set the hand to nan to make sure untouched pixels remain that value and not zero, which is a valid HAND height.\n",
    "for k,p in enumerate(tqdm(polygons)):\n",
    "    verbose=False    \n",
    "    mask,tf,win = rasterio.mask.raster_geometry_mask(dem, [p], crop=True, pad=True, pad_width=pad_width) #add 1 pixel. calculate_hand needs it.  \n",
    "    if win.width==1 or win.height==1: #padding may require this limit to change. \n",
    "        continue # The DEM is a thin line, skip this patch\n",
    "    not_mask=np.bitwise_not(mask)\n",
    "    #if polygon_ids[k] == 4120928640: #k=15 for polygon_ids, in hybas_id[70883]\n",
    "    #    verbose=True\n",
    "    if dem_nodata_mask[win.row_off:win.row_off+win.height,win.col_off:win.col_off+win.width].all():\n",
    "        continue # do not process if the entire polygon is nodata\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
    "        h=calculate_hand(np.squeeze(dem.read(window=win)), tf, pyproj.Proj(init=dem.crs.to_string()), mask=not_mask, verbose=verbose, acc_thresh=accumulation_threshold)\n",
    "    clip_hand=hand[win.row_off:win.row_off+win.height,win.col_off:win.col_off+win.width] #By reference\n",
    "    clip_hand[not_mask]=h[not_mask]\n",
    "hand[dem_nodata_mask]=nodata_fill_value\n",
    "print_duration(t)\n",
    "\n",
    "#this is for debugging only\n",
    "#hand_file=os.path.splitext(dem_file)[0]+f\"_hand_wNaN_{version.replace('.','_')}.tif\"\n",
    "#gdal_write(hand, dem_gT, filename=hand_file, srs_proj4=dem_proj4, nodata=np.nan, data_format=gdal.GDT_Float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b>3.4 Final Check for NaNs in HAND</b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\"> We check here again if there is any NaN values in HAND. This may be due to the fill_nan kernel being too small to fill all the NaNs, or alternatively if the basins do not cover the entire DEM due to errors in shapefile which may result in gaps.\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\"> Assuming we have a very large HAND that we can not fill-nan in one shot, we do this in smaller chunks. We find where the NaNs are, and separate non-touching NaNs into labelled objects using `ndimage.label()`. We then use this to loop over the HAND and fill NaNs over these small chunks. We again use assignment by reference to assign smaller nan-filled arrays into the large HAND array.\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\"> Finally, HAND can include NaNs due to coastal areas. We mask that using the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG). \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Create a Land Mask so that we do not try to NAN-Fill over the Ocean. \n",
    "t=time.time()\n",
    "if np.any(np.isnan(hand)):\n",
    "    print(f'{np.sum(np.isnan(hand))} NaN Pixels Detected in hand_result')\n",
    "    #generate nan_mask\n",
    "    #hand_type=hand.dtype\n",
    "    #hand_orig=hand.copy()\n",
    "    nan_mask=np.isnan(hand)\n",
    "    # Download GSHHG    \n",
    "    gshhg_zipfile=os.path.join(gshhg_dir,\"gshhg-shp-2.3.7.zip\")\n",
    "    gshhg_file=os.path.join(gshhg_dir,\"GSHHS_shp/f/GSHHS_f_L1.shp\")\n",
    "    if not os.path.exists(gshhg_zipfile) and not os.path.exists(gshhg_file):\n",
    "        #!wget -O {gshhg_zipfile} http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-shp-2.3.7.zip\n",
    "        urllib.request.urlretrieve(gshhg_url, gshhg_zipfile)\n",
    "    if not os.path.exists(gshhg_file): \n",
    "        #!unzip {gshhg_zipfile}\n",
    "        with zipfile.ZipFile(gshhg_zipfile, 'r') as zip_ref:          \n",
    "            zip_ref.extractall(path=gshhg_dir)\n",
    "    #if needed warp gshhg\n",
    "    if dem.crs.to_string() != 'EPSG:4326':\n",
    "        print(\"DEM and GSHHG projections differ.\")\n",
    "        #read extent shapes\n",
    "        gshhg_df=gpd.read_file(gshhg_file)\n",
    "        shapes=fiona_read_vectorfile(gshhg_file)\n",
    "        #find intersecting shapes\n",
    "        polygons=intersect(shapes, dem_poly_wgs84)              \n",
    "        gshhg=[]\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
    "            for p in polygons:\n",
    "                pt=transform_polygon(p.exterior.coords, s_srs='epsg:4326', t_srs='epsg:'+str(dem.crs.to_epsg()))\n",
    "                gshhg.append(point_coordinates_to_geometry(pt.exterior.coords))\n",
    "    else:\n",
    "        gshhg=fiona_read_vectorfile(gshhg_file)       \n",
    "    #generate land_mask for the DEM\n",
    " \n",
    "    land_mask,tf,win = rasterio.mask.raster_geometry_mask(dem, gshhg, crop=False, invert=True) #invert=If False (default), mask will be False inside shapes and True outside\n",
    "    #set ocean/sea values in hand to epsilon\n",
    "    if nodata_fill_value is not None:\n",
    "        hand[np.invert(land_mask)]=nodata_fill_value #sea_mask=np.invert(land_mask)\n",
    "    #find nan areas that are within land_mask    \n",
    "    joint_mask=np.bitwise_and(nan_mask,land_mask)\n",
    "    mask_labels, num_labels=ndimage.label(joint_mask)\n",
    "    print(f\"Number of NaN areas to fill: {num_labels}\")\n",
    "print_duration(t)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b>3.4.1 Nan-Fill Loop</b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\"> For some larger DEMs the kernel can crash during the loop. Therefore we take a break here to export some files as backup and delete unused variables. \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Run the NAN-fill operation.\n",
    "t=time.time()\n",
    "try: #If ran multiple times these variables might not exist. Hence they are in a try/except block.\n",
    "    del shapes, polygons, gshhg\n",
    "    del land_mask, nan_mask\n",
    "except:\n",
    "    pass\n",
    "print('Filling NaNs')\n",
    "\n",
    "#new nan_fill needs DEM. Might be better to NOT load it in the memory See: https://rasterio.readthedocs.io/en/latest/topics/windowed-rw.html\n",
    "demarray=dem.read(1)\n",
    "\n",
    "if np.any(np.isnan(hand)):\n",
    "    object_slices=ndimage.find_objects(mask_labels)    \n",
    "    tq=tqdm(range(1,num_labels))\n",
    "    for l in tq: #Skip first, largest label.        \n",
    "        #ids = np.argwhere(mask_labels==l)\n",
    "        #min0=max(ids[:,0].min()-1, 0)\n",
    "        #max0=min(ids[:,0].max()+1, mask_labels.shape[0])\n",
    "        #min1=max(ids[:,1].min()-1, 0)\n",
    "        #max1=min(ids[:,1].max()+1, mask_labels.shape[1])        \n",
    "        slices=object_slices[l-1] #osl label=1 is in object_slices[0]\n",
    "        min0=max(slices[0].start-1,0) \n",
    "        max0=min(slices[0].stop+1, mask_labels.shape[0])\n",
    "        min1=max(slices[1].start-1, 0)\n",
    "        max1=min(slices[1].stop+1, mask_labels.shape[1])       \n",
    "        mask_labels_clip=mask_labels[min0:max0, min1:max1]\n",
    "        h=hand[min0:max0, min1:max1] #by reference\n",
    "        d=demarray[min0:max0, min1:max1]\n",
    "        m=joint_mask[min0:max0, min1:max1].copy()\n",
    "        m[mask_labels_clip!=l]=0 #Maskout other flooded areas (labels) for this area. Use only one label.\n",
    "        if np.size(m)>1e6:\n",
    "            num_nan=m.sum()\n",
    "            tq.set_description(f\"Size: {num_nan}\")\n",
    "            if num_nan<1e6:\n",
    "                #hf=fill_nan(h.copy()) #break reference\n",
    "                hf=fill_nan_based_on_DEM(h.copy(), d.copy())\n",
    "                h[m]=hf[m] #copy nanfill by reference\n",
    "            else:\n",
    "                print(f'Filling {num_nan} pixels')\n",
    "                print('This can take a long time...')\n",
    "                hf=fill_nan_based_on_DEM(h.copy(), d.copy()) #break reference\n",
    "                h[m]=hf[m] #copy nanfill by reference\n",
    "        else:\n",
    "            hf=fill_nan_based_on_DEM(h.copy(), d.copy()) #break reference\n",
    "            h[m]=hf[m] #copy nanfill by reference\n",
    "print_duration(t)            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b><i>4. Export Raster </i></b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\">Output the HAND raster as a gdal generated geotiff. If a different format is preferred, can be provided to the gdal_write function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Write HAND to a file. \n",
    "hand_file=os.path.splitext(dem_file)[0]+f\"_hand_{version.replace('.','_')}.tif\"\n",
    "gdal_write(hand, dem_gT, filename=hand_file, srs_proj4=dem_proj4, nodata=np.nan, data_format=gdal.GDT_Float32)\n",
    "#cleaning up \n",
    "temporary_folder_object.cleanup()\n",
    "print('Overall process:')\n",
    "print_duration(t0)\n",
    "\n",
    "#Below is not working because file is outside Jupyter directory. See: https://stackoverflow.com/questions/24437661/retrieving-files-from-remote-ipython-notebook-server\n",
    "#local_file = FileLink(hand_file, result_html_prefix=\"Click here to download: \")\n",
    "#display(local_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"4\"> <b><i>4.1 Display HAND</i></b></font>\n",
    "<br><br>\n",
    "<font face=\"Calibri\" size=\"3\">Display the HAND output using std.dev. based min,max for colobar. Holes (no-data, not-a-number,nan) appear as white. We will not go below 0 since height above nearest drainage should not be negative below zero. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#Show HAND if plotting is requested. \n",
    "if show_plots:\n",
    "    #calculate mean and std.dev. for HAND\n",
    "    m=np.nanmean(hand) \n",
    "    s=np.nanstd(hand)\n",
    "    #minimum value should be no lower than 0. \n",
    "    hmin=max(m-2*s,0)\n",
    "    #If there is a large variation, just show the first 15m. \n",
    "    hmax=min(m+2.5*s,15)\n",
    "    print(f\"Setting colorbar limits as min:{hmin} max:{hmax}\")\n",
    "    pl.matshow(hand);pl.colorbar();pl.clim([hmin,hmax]);pl.title(f'Height Above Nearest Drainage {version} (m)');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font face=\"Calibri\" size=\"2\" color=\"gray\"> <i> Contacts: Batu Osmanoglu, MinJeong Jo\n",
    "    \n",
    "<b>Change Log</b> <br>\n",
    "2022/02/25: v0.1.13<br>\n",
    "-BugFix: Updating code to work with latest pysheds. <br>\n",
    "2021/08/18: v0.1.12<br>\n",
    "-BugFix: Improved interpolation with interpolating the DEM-HAND instead of HAND directly. We are referring to this as Height of Nearest Drainage (HOND=DEM-HAND). <br>\n",
    "2021/03/30: v0.1.11<br>\n",
    "-Feat: Implemented DEM nodata mask, to avoid calculation of HAND when the DEM has nodata values. These values are filled with nodata_fill_value (renamed ocean_fill_value).  <br>   \n",
    "2021/03/29: v0.1.10<br>\n",
    "-BugFix: Franz Meyer and batu adding ocean_fill_value to avoid long nan-fill operation in HYDRO30Workflow-v1. <br>\n",
    "2021/01/27: v0.1.9 <br>  \n",
    "-BugFix: Franz Meyer (UAF) reported an issue stating fiona_write_vectorfile() is failing in section 3.4 by giving a `fiona.errors.CRSError`. Implemented a workaround to avoid the function.<br> \n",
    "    2021/01/20: v0.1.8 <br>\n",
    "-Feat: Added zero-fill for low-lying nan areas in the HAND. These are often drainage areas with vegetation (or levees) around, denying any water from the basin to drain, despite being the lowest region in DEM. <br>\n",
    "2021/01/14: v0.1.7 <br>    \n",
    "-BugFix: Improved nan-handling for HAND calculation resolve_flats step. <br>\n",
    "-Feat: Added automatic accumulation threshold calculation. <br>\n",
    "2021/01/07: v0.1.6 <br>\n",
    "-Feat: Folder cleanup is added.<br>\n",
    "-Feat: Using zipfile and urllib instead of unzip and wget system calls. <br>\n",
    "2020/11/16: v0.1.5  <br>\n",
    "-BugFix: The EU HydroBasins link was set to AR<br>\n",
    "2020/10/16: 0.1.4 <br>\n",
    "-BugFix: Output was not respecting requested format, and was defaulting to 64bit. Also changed gdal_write() to generate Cloud-Optimized-Geotiff by default.  <br>\n",
    "-Feat: Performance upgrade to nan-filling. <br>\n",
    "2020/09/17: <br>\n",
    "-BugFix: Evan Smith (BYU) reported an issue which was resulting from the basin mapping to a 1pixel-wide DEM. Modified code to skip single pixel wide or high DEM patches. <br>\n",
    "2020/07/21: <br>\n",
    "-BugFix: Certain DEM projections was causing Inf values when global datasets were reprojected. Added a clipping step before projecting global datasets to DEM projection. <br>\n",
    "2020/06/02: <br>\n",
    "-BugFix: Added individual HydroBasin links for each region. <br>\n",
    "-BugFix: Clearing no-longer used variables before final nan-fill loop. Also skipping nan-fill if region is bigger than 1 Million pixels. <br>\n",
    "-BugFix: Moved export before display to make sure Hand is saved before attempting anything else when result is ready.<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hydrosar",
   "language": "python",
   "name": "conda-env-.local-hydrosar-py"
  },
  "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}