HAND and Flood Emergency Response

This Jupyter notebook illustrates the HAND workflow and its use in example flood emergency scenarios. The study area is Onion Creek (HUC10 code 1209020504).

This is also a demonstration of conducting geospatial anlysis with opensource toolkits (gdal) using an online Jupyter interface.

Environment required: CyberGIS-Jupyter for Water

  • We use CyberGIS' accelerated TauDEM version for d8 and d$\infty$ flow direction calculation.

    More about TauDEM

  • The USGS 3DEP Elevation dataset (a.k.a. National Elevation Dataset) is deployed on ROGER as a VRT raster.

  • The NHDPlus and associated water boundary dataset (WBD) are in Esri FileGDB format.

1. Data preparation and pre-processing

In [ ]:
# Program file path
import os
data_dir = '/home/jovyan/shared_data/data/onion-hand'
find_inlets = os.path.join(data_dir, 'find_inlets/build/find_inlets_mr')
    
# Format and performance parameters
np=os.cpu_count()
dsepsg="EPSG:4269"
bufferdist='0.2'

  Fetch the hands-on example data, from GIS in Water Resources (USU CEE6440, UT Austin CE 394K.3)

In [ ]:
!cp -r {data_dir} data

  The study area is Onion Creek in Texas. This watershed is a HUC10 unit with code 1209020504.

In [ ]:
!ogrinfo data/OnionHand.gdb
In [ ]:
!gdalinfo data/Onion3.tif

More about ogrinfo and gdalinfo

  To visualize a GeoTIFF, we need to:

  • Render color from values.
  • Make tiles (pyramids) from the colored image.
In [ ]:
# Colorize using hillshade
!gdaldem hillshade data/Onion3.tif Onion3-hillshade.tif
In [ ]:
# Make tiles with gdal2tiles (this is needed for visualization)
!gdal2tiles.py -t "Onion Creek DEM" -e -z 9-14 -a 0,0,0 -s epsg:4326\
    -r bilinear Onion3-hillshade.tif Onion3-hillshade

More about gdaldem and gdal2tiles.py

  To visualize a Shapefile, we only need to convert it to GeoJson

In [ ]:
# Extract flowline.shp from GeoDatabase (which will be used later)
!ogr2ogr -f "ESRI Shapefile" flowline.shp data/OnionHand.gdb NHDFlowline 2>/dev/null
In [ ]:
# Convert the Shapefile to GeoJson for visualization
!ogr2ogr -f geojson flowline.json flowline.shp

More about ogr2ogr

In [ ]:
from cybergis import Floret
Floret('Onion Creek DEM', 'input')\
.addTMSLayer('DEM','Onion3-hillshade')\
.addGeoJson('Flowline','flowline.json').display()

2. Workflow for computing HAND

2.1 Find inlets

  Use the dangle operation to find channel heads (inlets) of streams. They will guide the DEM-based stream calculation.

In [ ]:
! which find_inlets
In [ ]:
!find_inlets_mr -flow flowline.shp -dangle flow-inlets0.shp
!ogr2ogr -t_srs "EPSG:4269" flow-inlets.shp flow-inlets0.shp
In [ ]:
!ogr2ogr -f geojson flow-inlets.json flow-inlets.shp
In [ ]:
from cybergis import Floret
Floret('Flow Inlets', 'inlet').addTMSLayer('DEM','Onion3-hillshade')\
.addGeoJson('Flowline', 'flowline.json')\
.addGeoJson('Inlets','flow-inlets.json')\
.display()

  Rasterize inlet points to the same spatial extent as the input DEM. The output is the weight grid.

In [ ]:
ls data/getRasterInfo.py
In [ ]:
# Getting format info of the reference tif
import sys
names='fsizeDEM colsDEM rowsDEM nodataDEM xmin ymin xmax ymax cellsize_resx cellsize_resy'.split()
values=!{sys.executable} data/getRasterInfo.py data/Onion3.tif
print(zip(names, values))
for name,value in dict(zip(names,values)).items():
   para='%s="%s"'%(name,value)
   print (para)
   exec(para)
In [ ]:
!gdal_rasterize -ot Int16 -of GTiff -burn 1 -tr {cellsize_resx} {cellsize_resy} -te {xmin} {ymin} {xmax} {ymax} flow-inlets.shp weights.tif

More about gdal_rasterize

2.2 Pitremove

  Taudem pitremove: hydro-condition the DEM

In [ ]:
!mpirun -np {np} pitremove -z data/Onion3.tif -fel fel.tif

2.3 Compute flow direction

  Taudem d$\infty$: calc flow routing info that is used later to find the nearest stream from each cell

In [ ]:
!mpirun -np {np} dinfflowdir -fel fel.tif -ang ang.tif -slp slp.tif

 TauDEM d8: calc d8 flow directions for the computing of a stream grid weighted by the weight grid

In [ ]:
!mpirun -np {np} d8flowdir -fel fel.tif -p p.tif -sd8 sd8.tif

 TauDEM aread8: calc contributing area with the weight grid. i.e., only the contributing area starting from the channel heads denoted in the weight grid

In [ ]:
!mpirun -np {np} aread8 -p p.tif -ad8 ssa.tif -wg weights.tif -nc

 From the contributing area, extract streams areas (cellvalue = 1)

In [ ]:
!mpirun -np {np} threshold -ssa ssa.tif -src src.tif -thresh 1
In [ ]:
!gdaldem color-relief src.tif data/src.clr src-color.tif -of GTiff -alpha -nearest_color_entry
!gdal2tiles.py -e -z 9-14 -a 0,0,0 -s epsg:4326 -r bilinear src-color.tif src-color
In [ ]:
from cybergis import Floret
Floret('DEM-Derived streamline','src')\
.addTMSLayer('Derived Streamline', 'src-color')\
.addGeoJson('Flowline','flowline.json').display()

2.3 Calculate HAND

 With the d$\infty$ flowing routing and a stream grid that we know water flows on, calc HAND. We use vertical distance here; but horizontal distance or the combination can be used.

In [ ]:
!mpirun -np {np} dinfdistdown -fel fel.tif -ang ang.tif -src src.tif -dd hand.tif -m ave v
In [ ]:
# Colorization and tiling
!rm -rf ang-color
!gdaldem color-relief ang.tif data/ang.clr ang-color.tif -of GTiff -alpha
!gdal2tiles.py -e -z 9-14 -a 0,0,0 -s epsg:4326 -r bilinear ang-color.tif ang-color
In [ ]:
from cybergis import Floret
Floret('Dinf Distancedown','ang')\
.addTMSLayer('Dinf', 'ang-color')\
.addGeoJson('Flowline','flowline.json').display()
In [ ]:
# Colorization and tiling
!gdaldem color-relief hand.tif data/HAND-blues.clr hand-color.tif -of GTiff -alpha
!gdal2tiles.py -e -z 9-14 -a 0,0,0 -s epsg:4326 -r bilinear hand-color.tif hand-color
In [ ]:
from cybergis import Floret
Floret('Onion Creek HAND', 'onion-hand')\
.addTMSLayer('DEM', 'Onion3-hillshade')\
.addTMSLayer('HAND', 'hand-color')\
.addGeoJson('Flowline', 'flowline.json')\
.display()

3. Application

  Whose homes will be flooded if a 5m flood take place?

In [ ]:
# Filtering the HAND tif with threshold value
!gdal_calc.py -A hand.tif --outfile=depth.tif --calc="A<=5"
In [ ]:
# Convert the result into polygons
!gdal_polygonize.py depth.tif depth.shp # ed Becky

More about gdal_calc.py and gdal_polygonize.py

In [ ]:
!ogrinfo -al -so depth.shp
In [ ]:
# Filter the target areas from background
!ogr2ogr filtered_depth.shp depth.shp -where "DN=1"
In [ ]:
# Convert the result to a GeoJson
!ogr2ogr -f geojson filtered_depth.json filtered_depth.shp
In [ ]:
# The catchment of interest and the addressPoints in that catchment
!ogr2ogr -f geojson catchment.json data/catchment.shp
!ogr2ogr -f geojson addressPoints.json data/addressPoints.shp
In [ ]:
from cybergis import Floret
m = Floret('Catchment addresses and 5m-flood', '5m')
m.addTMSLayer('HAND', 'hand-color')
m.addGeoJson('Vunerable Areas', 'filtered_depth.json')
m.addGeoJson('Catchment', 'catchment.json')
m.addGeoJson('Address Points', 'addressPoints.json')
m.display()
In [ ]:
from cybergis import Editor
area = Editor('filtered_depth.shp')
addr = Editor('data/addressPoints.shp')
result=[]
for i in range(addr.lens):
    if -1 != area.index_of_first_feature_contains_point(*addr.shape(i).points[0]):
        result.append(i)
ans=addr.subshp(result)
ans.save('vunerableAddr')
!ogr2ogr -f geojson vunerableAddr.json vunerableAddr.shp
In [ ]:
from cybergis import Floret
m = Floret('Flooded addresses', 'flooded')
m.addTMSLayer('HAND', 'hand-color')
m.addGeoJson('Vunerable Areas', 'filtered_depth.json')
m.addGeoJson('Catchment', 'catchment.json')
m.addGeoJson('Vunerable Addresses', 'vunerableAddr.json')
m.display()

Congratulations on reaching the end!