NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Flow routing and watershed analysis

Resources:

Text file with outlet point coordinates:

Start GRASS GIS

Create a new mapset (called e.g. HW_hydrology) in nc_spm_08_grass7 project and change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter: 
cd
Download all text files with site locations (see above) to the selected directory. Now you can use the commands from the assignment requiring the text file without the need to specify the full path to the file. start GRASS

Compute flow direction, flow accumulation and subwatersheds

Compute flow direction, flow accumulation and subwatersheds with approx. size of 10000 cells from 30m NED. 
g.region raster=elev_ned_30m -p
r.watershed -s elevation=elev_ned_30m threshold=10000 accumulation=accum_10K drainage=draindir_10K basin=basin_10K

Extract more detailed streams from flow accumulation raster, and convert basins and streams to vector format: 

r.mapcalc "streams_der_30m = if(abs(accum_10K) > 100, 1, null())"
r.to.vect -s basin_10K output=basin_10K type=area
r.thin streams_der_30m output=streams_der_30m_t
r.to.vect -s streams_der_30m_t out=streams_der_30m type=line

Generate relief shaded map of basins and display extracted streams along with the official Wake county streams (red).
How do the derived streams compare with the official stream map? How can you modify the mapcalc expression to make stream origins fit better with the official stream map? 

d.erase
d.his hue=basin_10K intensity=elevation_shade brighten=40
d.vect basin_10K type=boundary
d.rast lakes
d.vect streams_der_30m color=blue
d.out.file mystreams
d.vect streams color=red
d.out.file streams_compare

Note that you can also use a MFD routing for r.watershed combined with r.stream.extract for vectorized streams with topology of stream network. Read the manual pages to learn more.

Create map of DEM depressions

Depression filling is sometimes necessary for certain flow routing algorithms but it can significantly alter the elevation data leading to srtificial straight channels. Find out how extensive the depressions are in our DEM.
Note that r.watershed tool doesn't need any depression filling thanks to its underlying algorithm which uses least cost path to get over depressions. 
g.region raster=elevation -p
r.fill.dir input=elevation output=elev_fill1 direction=dir1 areas=unres1
r.fill.dir input=elev_fill1 output=elev_fill2 direction=dir2 areas=unres2
r.fill.dir input=elev_fill2 output=elev_fill3 direction=dir3 areas=unres3
r.mapcalc "depr_bin = if((elevation-elev_fill3) < 0., 1, null())"
r.colors depr_bin color=blues

Remove all previously used layers from the Layer Manager and display the new results, compare the derived depressions with actual lakes: 

d.erase
d.rast elevation
d.vect roadsmajor
d.rast depr_bin
d.vect lakes type=area fill_color=blue
d.out.file depressions
What is the area mapped as depression in [ha] and [%] ? (hint: use r.report with depr_bin)
Optional: How would you compute the total volume of mapped depressions?

Derive contributing area for a given outlet

Set region to the high resolution study area and zoom to it, compute flow accumulation, flow direction and basins with 15000 cells threshold: 
g.region rural_1m -p
r.watershed -as elevation=elev_lid792_1m threshold=15000 accumulation=accum_15K drainage=draindir_15K basin=basin_15K

Remove previous layers. Display extracted streams over aerial image: 

d.erase
d.rast ortho_2001_t792_1m
d.rast accum_15K values=1000-1000000

Delineate contributing area (watershed) associated with an outlet located on the extracted stream and convert it to vector format. Import coordinates of the outlet to include it in the map in the next step. 

r.water.outlet input=draindir_15K output=basin_A30 coordinates=638845.52,220085.26
r.to.vect -s basin_A30 output=basin_A30 type=area
v.in.ascii input=outlet_point.txt output=outletA30 separator=space

Display watershed boundary along with contours and compute the watershed area (category 1)

d.erase
d.rast ortho_2001_t792_1m
d.vect basin_A30 type=boundary color=green width=2
r.contour elev_lid792_1m output=elev_lid792_cont_1m step=1 minlevel=104
d.vect elev_lid792_cont_1m color=white
d.vect outletA30 color=yellow size=15 icon=basic/circle
d.out.file watershedA30
r.report basin_A30 unit=h,a
What is the delineated watershed area [ha] ?
Optional: How would you compute average slope of this contributing area?

Compute refined flow pattern using D-inf

Compare upslope and downslope flow lines: On what type of landform (ridge, valley) do the upslope and dowslope flowlines converge? 
g.region raster=elev_lid792_1m -p
r.flow elevation=elev_lid792_1m flowline=flowlines flowlength=flowlg_1m flowaccumulation=flowacc_1m
r.flow -u elevation=elev_lid792_1m flowlength=flowlgup_1m flowaccumulation=flowaccup_1m

Display maps along with contours to see the relation of flow pattern to terrain: 

d.erase
d.rast flowacc_1m
d.vect elev_lid792_cont_1m color=red
d.vect flowlines
d.out.file flowdown
d.erase
d.rast flowaccup_1m
d.vect elev_lid792_cont_1m color=red
d.out.file flowlines

Compare the multiple flow direction flow accumulation with the D8 and Dinf

Compute flowaccumulation using the default MFD method and display the flowaccumulation maps from MFD and SFD (the Dinf maps flowacc_1m was displayed above): 
r.watershed -a elevation=elev_lid792_1m threshold=15000 accumulation=accum_mfd_15K
d.rast accum_mfd_15K
d.vect streams co=red
d.out.file stream_mfd
d.erase
d.rast accum_15K
d.vect streams co=red
d.out.file stream_sfd
Compare the MFD and SFD results of r.watershed. Why are they different, which is more realistic and why?

Compute wetness index

g.region rural_1m -p
r.topidx elev_lid792_1m output=wetness_1m
r.colors map=wetness_1m color=water
d.rast wetness_1m
d.out.file wetness
Optional: Compute the wettness index using the formula in the lecture (hint: you will need to compute slope using elev_lid792_1m and use accum_mfd_15K as contributing area. Check the values by comparing your output with wetness_1m computed with r.topidx.

Display elev_lid792_1m in 3D and drape over wetness_1m as color.
Note: switch off all layers except for elev_lid792_1m.

Create a map of flooded area

Create a map of flooded area for a given water level and seed point: 
g.region rast=elev_lid792_1m -p
r.lake elevation=elev_lid792_1m water_level=113.5 lake=flood1 coordinates=638728,220278
d.rast elev_lid792_1m
d.rast flood1
d.out.file floodedarea

Increase water level to 113.7m and 115.0m and create flooded area maps at these two levels. Include the maps of flooded areas at these water levels in your report. Why is the flooded area much larger than for 113.5m?

Optional: Assess and mitigate impact of the road on flowrouting.

Compare the extracted streams (accum > 1500) with official stream data: 

g.region rast=elev_lid792_1m -p
d.rast ortho_2001_t792_1m
d.rast accum_5K values=1500-1000000
d.vect streets_wake color=red
d.vect streams color=green
d.out.file streamcompare

Carve a pre-defined channel given by the stream data into DEM

r.carve raster=elev_lid792_1m vector=streams width=2 depth=0.8 output=elev_lidcarved_1m
r.colors elev_lidcarved_1m raster=elev_lid792_1m
d.rast elev_lidcarved_1m

Extract streams from the carved DEM and compare with the official streams map.
What is the difference between accum_5K, accumc_5K1m and streams?
Explain the advantage and disadvantage of carving. 

r.watershed -as elevation=elev_lidcarved_1m accumulation=accumc_5K1m
d.rast accumc_5K1m values=1500-10000000
d.vect streams
d.out.file streamcarved