GIS/MEA582: Geospatial Modeling and Analysis

Flow routing and watershed analysis

Resources:

GRASS GIS overview and manual
Recommendations and tutorial how to use GUI from the first assignment

Text file with outlet point coordinates:

point coordinates (outlet_point.txt)

Start GRASS GIS

In startup pannel set GRASS GIS Database Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata. For GRASS Location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for GRASS Mapset create a new mapset (called e.g. HW_hydrology) and click Start GRASS session.

grass83

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:

r.mapcalc "streams_der_30m = if(abs(accum_10K) > 100, 1, null())"

Convert to vector format and display results along with the official Wake county streams (red):

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 shaded map and display: 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.vect streams color=red
d.out.file mystreams

Create map of DEM depressions

Depression filling is often necessary for certain flow routing algorithms but it can alter the elevation data significantly. Find out how extensive the depressions are in our DEM.
Note that r.watershed 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

Derive contributing area for a given outlet

Set region to the high resolution study area and zoom to it:

g.region rural_1m -p
r.watershed -as elevation=elev_lid792_1m threshold=5000 accumulation=accum_5K drainage=draindir_5K basin=basin_5K

Remove previous layers. Display extracted streams over aerial image:

d.erase
d.rast ortho_2001_t792_1m
d.rast accum_5K values=1500-1000000

Identify outlet on the extracted stream.

Create a vector map with the point east=638845.52 north=220085.26 (download the text file) that has accum_5K=224510.

v.in.ascii input=outlet_point.txt output=outletA30 separator=space
d.vect outletA30 color=yellow

Delineate the contributing area associated with this outlet and convert it to vector format:

r.water.outlet input=draindir_5K output=basin_A30 coordinates=638845.52,220085.26
r.to.vect -s basin_A30 output=basin_A30 type=area

Display watershed boundary along with contours.
Compute the watershed area (category 1)

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.out.file watershedA30
r.report basin_A30 unit=h,a

Assess and mitigate impact of the road on flowrouting

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

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

Compute refined flow pattern using D-inf

Compare upslope and downslope flow lines: On what type of landform (ridge, valley) they 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 relation 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 result with D8 and Dinf

Compare the results of r.watershed: accum_5K and accum_mfd_5K1m. Why are they different, which is more accurate and why?

r.watershed -a elevation=elev_lid792_1m threshold=5000 accumulation=accum_mfd_5K1m
d.rast accum_mfd_5K1m
d.vect streams
d.out.file stream_mfd
d.erase
d.rast accum_5K
d.vect streams
d.out.file stream_sfd

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

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_5K1m 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:

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.