4

I am new to spatial analysis and even newer to hydrology. I have a 30arc DEM elevation raster and I was wondering if I can use it as an input to get the approximate presence of rivers in a geographically determined zone. I know I can use rivers shapefile to do this, but since I am working with cells as my unit observation, it is cleaner get the data from the DEM raster.

raster package in R seems appropriate to this. Nonetheless, it only calculates hill, slope and roughness.

5
#You can run GRASS 6.4 commands in R using the spgrass6 package    
require(spgrass6)

# connect to existing GRASS layer location using Mac OS X
loc <- initGRASS("/Applications/GRASS-6.4.app/Contents/MacOS",
                 gisDbase="/Users/junf/Documents/WatershedModel/MyDataFolder", 
                 location="colonia", mapset="flood",override=T,
                 SG="fil_srtm30_fill_dist2out")

# use Grass 6.4 r.watershed to delineate streams
execGRASS(cmd='r.watershed', flags='overwrite', 
    parameters =  list(elevation='val_srtm30@flood', threshold='2000', 
          drainage='wat30_drain', accumulation='wat30_accu', basin='wat30__basin', 
          stream='wat30__stream', memory='2000'))
# thin streams to make them 1-pixel wide
execGRASS('r.thin', parameters = list(input='wat30__stream', output='wat30__stream')) 
# convert stream raster to vector
execGRASS('r.to.vect', parameters = list(input='wat30__stream', output='v_stream_wat30'))

#The final output stream vector is v_stream_wat30 which remains GRASS format
3
  • There is a spgrass7 library as well. – Jeffrey Evans Feb 26 '16 at 4:51
  • Jun, thanks for your answer. Can you be more explicit in your response? What are those lines of code doing? I understand the parameters in the spgrass package, but which files do I need to get the watershed or rivers flow? – topcat Feb 29 '16 at 17:09
  • @Iván Higuera Mendieta please check my answer below with a reproductible code and sample data – Jun Folledo Mar 1 '16 at 8:15
4

Two free software solutions for watershed and stream delineation are the GUI software SAGA GIS (which can also be accessed through the sextante toolbox in QGIS) and David Tarboton's TauDEM. You can also perform comprehensive hydrological modeling, watershed and stream delineation using GRASS GIS but it is a much steeper learning curve. That said, GRASS tools can be accessed directly in QGIS GUI, through the plugins, and provides a very nice interface to GRASS.

Both SAGA and TauDEM are very good solutions, and I would say much better and more flexible than ArcGIS. SAGA provides a variety of DEM conditioning options as well as multiple algorithms for deriving flow direction. The TauDEM software is specifically designed for hydrologic modeling and stream delineation with a clear workflow.

One note on deriving the required rasters for stream delineation (flow direction and flow accumlation). The traditional D8 flow direction algorithm often results in unrealistic flow routes (square channels) and in my experience the Dinf (D infinity) flow direction algorithm yields superior results to D8. The flow direction raster will be used to derive flow accumulation which is then the raster utilized for channel delineation.

1

I am adding a reproductible code, including a GRASS database (testR) so you can run it as long as you have installed GRASS 6.4 and the spgrass6 package in R

The GRASS database which contains 1 DEM raster file named strm30_utm. Copy the testR folder to any folder, in my code, I copied it at the Documents "/Users/junf/Documents"

GRASS database: https://drive.google.com/file/d/0B5upqlYhcHH4ZzdkQkVoX1Mtbmc/view?usp=sharing

require(spgrass6)

# create connection to GRASS database
initGRASS("/Applications/GRASS-6.4.app/Contents/MacOS", gisDbase="/Users/junf/Documents/", location="testR", mapset="wtshd",override=T)

unlink_.gislock() # unlock to be able to open testR database in GRASS

# create subwatershed: 
# input DEM is strm30_utm (stored in attached Grass database), 30 meter resolution
# output stream raster is up_stream
execGRASS('r.watershed', flags='overwrite',parameters = list(elevation='strm30_utm', threshold=2000, stream='up_stream', basin='r_basin'))

# r.thin thins non-NULL pixels so that each line is only 1 pixel wide, required before converting to vector
execGRASS('r.thin', flags='overwrite',parameters =  list(input='up_stream', output='r_strm_thin'))

# convert raster (r_strm_thin) to vector (v_stream)
execGRASS('r.to.vect', flags='overwrite',parameters = list(input='r_strm_thin', output='v_stream'))

# clean so that each line segment is at least 65 meters long 
execGRASS('v.clean', flags=c('overwrite'),parameters =list(input='v_stream', output='v_stream_clean', type='line',tool='snap',thresh=65))

# display output vector
system("d.mon x0")
system("d.rast r_basin")
system("d.vect v_stream_clean")

# export stream vector to shapefile
execGRASS('v.out.ogr', parameters=list(input='v_stream',type="line", dsn="/Users/junf/Documents/testR/",format="ESRI_Shapefile"))
1
  • Please format your code as code. Otherwise it's rather unreadable... – dru87 Mar 1 '16 at 7:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.