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¶
# Program file path
import os
data_dir = '/home/jovyan/shared_data/data/onion-hand'
# Format and performance parameters
np=4 #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)
!cp -r {data_dir} data
The study area is Onion Creek in Texas. This watershed is a HUC10 unit with code 1209020504.
!ogrinfo data/OnionHand.gdb
!gdalinfo data/Onion3.tif
To visualize a GeoTIFF, we need to:
- Render color from values.
- Make tiles (pyramids) from the colored image.
# Colorize using hillshade
!gdaldem hillshade data/Onion3.tif Onion3-hillshade.tif
# 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
# Extract flowline.shp from GeoDatabase (which will be used later)
!ogr2ogr -f "ESRI Shapefile" flowline.shp data/OnionHand.gdb NHDFlowline 2>/dev/null
# Convert the Shapefile to GeoJson for visualization
!ogr2ogr -f geojson flowline.json flowline.shp
More about ogr2ogr
from cybergis import Floret
Floret('Onion Creek DEM', 'input')\
.addTMSLayer('DEM','Onion3-hillshade')\
.addGeoJson('Flowline','flowline.json').display()
!find_inlets_mr -flow flowline.shp -dangle flow-inlets0.shp
!ogr2ogr -t_srs "EPSG:4269" flow-inlets.shp flow-inlets0.shp
!ogr2ogr -f geojson flow-inlets.json flow-inlets.shp
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.
ls data/getRasterInfo.py
# 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)
!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
!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
!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
!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
!mpirun -np {np} aread8 -p p.tif -ad8 ssa.tif -wg weights.tif -nc
From the contributing area, extract streams areas (cellvalue = 1)
!mpirun -np {np} threshold -ssa ssa.tif -src src.tif -thresh 1
!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
from cybergis import Floret
Floret('DEM-Derived streamline','src')\
.addTMSLayer('Derived Streamline', 'src-color')\
.addGeoJson('Flowline','flowline.json').display()
2.4 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.
!mpirun -np {np} dinfdistdown -fel fel.tif -ang ang.tif -src src.tif -dd hand.tif -m ave v
# 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
from cybergis import Floret
Floret('Dinf Distancedown','ang')\
.addTMSLayer('Dinf', 'ang-color')\
.addGeoJson('Flowline','flowline.json').display()
# 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
from cybergis import Floret
Floret('Onion Creek HAND', 'onion-hand')\
.addTMSLayer('DEM', 'Onion3-hillshade')\
.addTMSLayer('HAND', 'hand-color')\
.addGeoJson('Flowline', 'flowline.json')\
.display()
# Filtering the HAND tif with threshold value
!gdal_calc.py -A hand.tif --outfile=depth.tif --calc="A<=5"
# Convert the result into polygons
!gdal_polygonize.py depth.tif depth.shp # ed Becky
More about gdal_calc.py and gdal_polygonize.py
!ogrinfo -al -so depth.shp
# Filter the target areas from background
!ogr2ogr filtered_depth.shp depth.shp -where "DN=1"
# Convert the result to a GeoJson
!ogr2ogr -f geojson filtered_depth.json filtered_depth.shp
# 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
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()
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
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()