Watershed Analysis using GRASS GIS

After following an exercise provided at <http://www.ing.unitn.it/~grass/docs/tutorial_64_en/htdocs/esercitazione/dtm/dtm4.html>, I have been playing around with different thresholds and flow routing algorithms in GRASS GIS to try to reproduce the National Hydrography Dataset Plus provided by <http://www.horizon-systems.com/nhdplus/>. Below is the NHDPlus elevation, flowlines, and waterbodies data for Puerto Rico. I also added the HUC10 boundary for the Rio Grande de Arecibo basin.

Rio Grande de Arecibo basin in Puerto Rico with flowlines and waterbodies from NHDPlus

One thing I've been concerned with is representing the waterbodies accurately. From the GRASS GIS book, "The r.watershed module does not require filling of depressions (pits, sinks) in DEM prior to its application because it uses the least-cost algorithm to traverse the elevation surface to the outlet." So, I wanted to test some things.

Close up of the southwest corner of the watershed where Lago Garza is located

Two lakes in the eastern-central part of the watershed

Firstly, I ran r.watershed using SFD and MFD without filling in the depressions.

Using r.watershed with threshold of 10000 and single flow direction from the command line...

r.watershed elevation=elev_cm accumulation=rwater_accum10Ksfd   drainage=rwater_draindir10Ksfd basin=rwater_basin10Ksfd stream=rwater_streams10Ksfd threshold=10000

And, the results...

Basin map from using threshold of 10000 and single flow direction

Accumulation map with threshold of 10000 and single flow direction

Here is a close up of the southwest corner of the watershed

Close up of southwest corner accumulation map for the watershed

Taking the log of the accumulation map

r.mapcalc 'log_rwater_accum10Ksfd=log(abs(rwater_accum10Ksfd)+1)'

log accumulation

Here is the close up from using the MFD method

Accumulation map using MFD method

And taking the log again

Log accumulation map using MFD method

Here are the two lakes from the middle of the watershed

Log accumulation map of the two lakes in the middle of the watershed using MFD method

Next, I tried r.watershed after filling in the depressions.

r.fill.dir elev_cm el=elev_fill1 dir=dir1 areas=unres1

r.fill.dir elev_fill1 el=elev_fill2 dir=dir2 areas=unres2

r.fill.dir elev_fill2 el=elev_fill3 dir=dir3 areas=unres3

r.mapcalc 'depr_bin=if((elev_cm - elev_fill3)<0., 1, null())'

NHD Plus waterbodies outlined in blue and the depression filled raster map shown in fuschia

The accumulation map after running r.watershed with the depressions filled seems like it would represent the streams better.

Accumulation map when using depression map and MFD method

However, the resulting streams are not as expected. I tried using each level of filled depressions, and different ways to extract the streams from the accumulation map, but none gave good results.

Streams resulting from using a depression filled map

Here is a map with the streams derived from using r.watershed with SFD on top of the NHDPlus flowlines.

Streams derived from r.watershed using SFD in fuschia and the NHDPlus flowlines in blue

Here is SFD, MFD, and the NHDPlus all in one.

NHDPlus(blue), SFD(magenta), and MFD(light brown)

The MFD method using

r.mapcalc 'inf_rivers=if(log_accumulation>6)'

to extract the streams seems to match the NHDPlus flowlines the best.