Wednesday, June 29, 2011

Using r.buffer for watershed analysis

A problem that I ran into when using r.watershed with a masked basin was that the streams were being cut off, and the resulting basin was incorrect.

Close up of the HUC 10 watershed boundary in black and the streams derived from r.watershed in blue and pink


basin resulting from original mask


To solve this I created a buffer raster of the original raster being used as a mask.

r.buffer input=rioGrandeHUC10mask output=rioGrandeBuffer distances=100

Then, following the grass book

r.watershed -f elev_cm_rioGrande2 thresh=10000 accum=accum_10K drain=draindir_10K basin=basin_10K stream=rivers_10K --o
 
r.mapcalc "streams_der_30m=if(abs(accum_10K)>100,1,null())"

r.thin streams_der_30m out=streams_der_30m_tmp

r.to.vect -s streams_der_30m_tmp out=streams_der_30m

r.to.vect -s basin_10K out=basin_10K feature=area

Here is the resulting stream

Streams in green after using a buffer with 100m distance


r.water.outlet drainage=draindir_10K basin=outlet_O8 easting=170670 northing=270651 --o

r.mapcalc 'outlet_O8=if(outlet_O8==1,1,null())'


Resulting basin with buffered mask

Wednesday, June 22, 2011

More NHDPlus with GRASS GIS

Here is the flow accumulation grid from the NHDPlus dataset


Here, I take the log of the accumulation again:

r.mapcalc 'log_fac=log(abs(fac)+1)


Extract the streams and convert them to vector:

r.mapcalc 'inf_rivers=if(log_fac>6)'
r.thin input=inf_rivers out=thin_rivers
r.to.vect -s thin_rivers out=inf_rivers

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.

 

Sunday, June 19, 2011

Learning JGrassTools

I've been spending some time learning how to use JGrassTools. This is one of the examples from the wiki <http://code.google.com/p/jgrasstools/wiki/JGrassTools4CLI> that I started with. Hopefully, more posts will follow.


// creation of the simulation object
sim = new oms3.SimBuilder(logging:'ALL').sim(name:'pit+viewer') {
   // the model space
   model {
      // space for the definition of the acting components
      components  {
      'pitfiller'    'pit'
      'reader'        'rasterreader'
      'writer'        'rasterwriter'
      'viewer'        'mapsviewer'
      }
     
      // initialization of the parameters
      parameter  {
      'reader.file'        '/Users/justinjent/GrassData/puertoRico_GRS80_NAD83_UTM19/jgrass/cell/rioGrande'
      'writer.file'        '/Users/justinjent/GrassData/puertoRico_GRS80_NAD83_UTM19/jgrass/cell/rioGrande_pit'
      }
    
      // connection of the different components
      connect  {
      'reader.outRaster'    'pitfiller.inElev'
      'pitfiller.outPit'    'writer.inRaster'
      'pitfiller.outPit'    'viewer.inRaster'
      }
   }
}
// start of the simulation
sim.run();