generate-hand-notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Generate your own GLO-30 HAND tile\n",
    "\n",
    "This notebook will recreate a GLO-30 HAND tile, or let you generate a new one with using a different accumulation value.\n",
    "\n",
    "## 0. Install dependencies\n",
    "\n",
    "To run this notebook, you'll need the [ASF Tools](https://hyp3-docs.asf.alaska.edu/tools/asf_tools/) Python package, [numpy](https://numpy.org/), [rasterio](https://rasterio.readthedocs.io/), and [shapely](https://shapely.readthedocs.io/en/stable/manual.html). You can install them into a conda environment with:"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "%conda install -c conda-forge asf_tools numpy rasterio shapely"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 1. Generate HAND\n",
    "First, select a DEM tile and then provide the [virtual path](https://gdal.org/user/virtual_file_systems.html) to it"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "dem_tile = '/vsicurl/https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N41_00_W088_00_DEM/Copernicus_DSM_COG_10_N41_00_W088_00_DEM.tif'"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Extract some metadata from the DEM tile"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import rasterio\n",
    "\n",
    "with rasterio.open(dem_tile) as dem:\n",
    "        dem_bounds = dem.bounds\n",
    "        dem_meta = dem.meta"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Buffer DEM by 0.5 degrees to avoid clipping hydrobasin boundaries"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from shapely.geometry import GeometryCollection, box\n",
    "\n",
    "dem_geometry = box(*dem_bounds)\n",
    "dem_buffered = GeometryCollection([dem_geometry.buffer(0.5)])"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Build a temporary VRT of all DEM tiles that intersect our buffered geometry"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "from asf_tools.dem import prepare_dem_vrt\n",
    "\n",
    "buffered_dem_vrt = Path(dem_tile).name.replace('.tif', '_buffered.vrt')\n",
    "prepare_dem_vrt(buffered_dem_vrt, dem_buffered)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Calculate an initial HAND over the entire buffered area using the specified accumulation threshold"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from asf_tools.hand.calculate import calculate_hand_for_basins\n",
    "\n",
    "accumulation_threshold = 100\n",
    "\n",
    "hand_raster_buffered = Path(dem_tile).name.replace('DEM.tif', f'HAND_ACC{accumulation_threshold}_buffered.tif')\n",
    "calculate_hand_for_basins(hand_raster_buffered, dem_buffered, buffered_dem_vrt, acc_thresh=accumulation_threshold)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Crop out buffered HAND to the original DEM tile size"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "with rasterio.open(hand_raster_buffered) as sds:\n",
    "    window = rasterio.windows.from_bounds(*dem_bounds, sds.transform)\n",
    "    out_pixels = sds.read(\n",
    "        1, window=window, out_shape=(dem_meta['height'], dem_meta['width']),\n",
    "        resampling=rasterio.enums.Resampling.bilinear\n",
    "    )"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Mask the ocean pixels as identified in the WBM auxiliary DEM files"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "wmb_tile = dem_tile.replace('DEM/Copernicus', 'DEM/AUXFILES/Copernicus')\n",
    "wmb_tile = wmb_tile.replace('_DEM.tif', f'_WBM.tif')\n",
    "\n",
    "with rasterio.open(wmb_tile) as wbm:\n",
    "    wbm_pixels = wbm.read(1)\n",
    "\n",
    "out_pixels = np.ma.masked_where(wbm_pixels == 1, out_pixels)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Write out the final HAND tile"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from asf_tools.composite import write_cog\n",
    "\n",
    "hand_raster = Path(dem_tile).name.replace('DEM.tif', f'HAND_ACC{accumulation_threshold}.tif')\n",
    "write_cog(hand_raster, out_pixels, transform=dem.transform.to_gdal(), epsg_code=dem.crs.to_epsg())"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Check out the HAND you made!"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import rasterio.plot\n",
    "\n",
    "with rasterio.open(hand_raster) as hand:\n",
    "    rasterio.plot.show(hand, title=hand_raster, vmin=0, vmax=10, cmap='binary')"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}