Creating a Spatial Stream Network with GRASS GIS. Part I

17 Mar 2015

Recently, I wanted to use the SSN R package to predict along a stream network. To do this you need to set up the GIS files using the ArcGIS plugins here. I don’t have access to ArcGIS (and I don’t want it) because I only have access to ubuntu and OS X. This is why I decided to develop several GRASS GIS addons to provide the capability to create ssn objects to be used in R. This post will outline the first few steps used in delineating catchments in R and I will follow this post up with the addons later on.

I am going to outline how to create a ssn object in GRASS using the workshop data which can be found here: ssn-short-course.

Step 1: Creating the project

1. Create a new Project location
2. Give it a name
3. You now need to set the project, resolution and bounds of your project. The easiest way to do this is to select the third option and choose the dem file datasets/gisdata/dem/w001001.adf.
4. You will now be given the option to import the file into the project. Select No - will we do this later
5. Now you will be given the option to create a working directory. From my understanding this option exists to provide flexibility for multiple users. I prefer to use the term “working”, but this doesn’t matter.
6. Now you can start the grass session.

Step 2: Importing the maps

Below I will provide the list of commands which can be used directly in the command console within GRASS. You don’t need to use the commands as the GUI does have list most of the commands I will use. You can also use the command console to open an additional window for each grass command. For example you can type: r.in.gdal and press enter. This will open the import window.

Import the DEM map (the ‘e’ flag tells GRASS to set the bounds and resolution to the DEM map)

bash r.in.gdal input=/Users/jasonlessels/work/Projects/ssn-short-course/datasets/gisdata/dem/w001001.adf output=dem -e

You can display the DEM map

To import the developed land classified map

bash r.in.gdal input=/Users/jasonlessels/work/Projects/ssn-short-course/datasets/gisdata/developd/w001001.adf output=developed

Import the site location vector layer

bash v.in.ogr input=/Users/jasonlessels/work/Projects/ssn-short-course/datasets/gisdata/sub_sites.shp layer=sub_sites output=sites

Import the streams

bash v.in.ogr input=/Users/jasonlessels/work/Projects/ssn-short-course/datasets/gisdata/streams.shp layer=streams output=streams

Step 3: Filling in sinks and burning in streams.

Now I will outline how you can correct for sinks in the DEM raster, and if you want to, burn in the streams. Burning in the streams is handy if you have access to data with the locations of real streams.

Convert the streams (vector) to a raster

bash v.to.rast --overwrite input=streams@working type=line output=gov_streams use=cat

Subtract 10 m elevation from the DEM where streams lie.

bash r.mapcalc expression='dem_burn = if(!isnull(gov_streams), dem-10, dem)'

Fill the DEM and create a direction map.

bash r.fill.dir input=dem_burn@working output=dem_filled direction=fill_direction

Create accumulation map

bash r.watershed elevation=dem_filled threshold=250000 accumulation=accum

Extract streams using the r.streams modules.

bash r.stream.extract elevation=dem_filled accumulation=accum threshold=2500 stream_rast=stream_rast stream_vect=stream_vect direction=direction

Install the additional r.stream* addons.

bash g.extension extension=r.stream.order operation=add g.extension extension=r.stream.snap operation=add g.extension extension=r.stream.distance operation=add

Obtain the stream order

bash r.stream.order stream_rast=stream_rast direction=direction elevation=dem_filled accumulation=accum stream_vect=stream_order

The r.stream.order function add several additional columns with data. The names get truncated when they are exported as a shape file. The next few lines truncates the names now.

bash v.db.renamecolumn map=stream_order column='next_stream, nxt_str' v.db.renamecolumn map=stream_order column='scheidegger, scheide' v.db.renamecolumn map=stream_order column='source_elev, src_elev' v.db.renamecolumn map=stream_order column='outlet_elev, olt_elev' v.db.renamecolumn map=stream_order column='stream, rid' v.db.renamecolumn map=stream_order column='out_dist, upDist'

Snap the points to the streams using the r.stream.snap function. This will create a new vector layer. The radius is set to 10 m and the function will move the sites to the streams if the points are within 10 meters of a stream. This function also provides the ability to see if the points were moved or not.

bash r.stream.snap input=sites output=snap_sites stream_rast=stream_rast accumulation=accum radius=10

Add an attribute table to the snap_sites vector (Seems to be a bug)

bash v.db.addtable map=snap_sites table=snap_sites

Join the table with the sites table - this can take a minute or so.

bash v.db.join map=snap_sites column=cat otable=sites ocolumn=cat

Step 4. Site distance ratio

Now estimate the distance ratio for each. This is done by getting the distance to the end of each node for every cell and dividing it by the length of each edge.

Get the distance to the end of the entire network.

bash r.stream.distance -o stream_rast=stream_rast direction=direction elevation=dem_filled method=downstream distance=distance

Get the length of each edge and convert it to a raster.

bash v.to.rast input=stream_order type=line output=edge_length use=attr attrcolumn=length v.to.rast input=stream_order type=line output=upDist use=attr attrcolumn=upDist

Calculate the distance ratio for each edge.

bash r.mapcalc expression='distance_ratio=if(!isnull( edge_length),(distance - (upDist-edge_length))/edge_length , null())'

Calculate point upDist - this should be the same as upDist, but do it again anyway.

bash r.mapcalc expression='point_dist=if(!isnull( edge_length),(distance_ratio*edge_length)+(upDist-edge_length) , null())'

Step 5. Update the sites vector

Firstly, new columns have to be added to the snap_sites vector

bash v.db.addcolumn map=snap_sites columns='rid INT, netID INT, pid INT, locID INT, ratio double precision, NEAR_X double precision, NEAR_Y double precision, upDist double precision, elev double precision'

Now use the v.what.vect function to get the corresponding values from the raster maps.

bash v.what.vect map=snap_sites column=rid qmap=stream_order qcolumn=rid v.what.vect map=snap_sites column=netID qmap=stream_order qcolumn=netID v.what.rast map=snap_sites raster=distance_ratio column=ratio v.what.rast map=snap_sites raster=point_dist column=upDist v.what.rast map=snap_sites raster=dem column=elev

Create a prediction id column based on the

bash v.db.update map=snap_sites layer=1 column=pid qcolumn=cat v.db.update map=snap_sites layer=1 column=locID qcolumn=cat

Add x and y coordinates to the snap_sites vector using the v.to.db function.

bash v.to.db map=snap_sites option=coor col=NEAR_X,NEAR_Y

Step 6. Create a prediction grid

Create the predictions points using the v.to.points function. This will create a point every 1500 m along the stream

bash v.to.points input=stream_order type=line output=pred_pts dmax=1500

The vector created above has two tables - the table of stream_order and a new table for all the new points, named;

1. table 1 = pred_pts_1
2. table 2 = pred_pts_2

The next step is to use the lcat (which is the same as the rid) to merge the tables together.

bash v.db.join map=pred_pts layer=2 column=lcat otable=pred_pts_1 ocolumn=rid scolumns=rid,length,upDist,netID

Calculate the up stream distance for each point using the along and upDist columns. The v.to.points function added an along column, which is the distance from the node of every edge.

bash v.db.update pred_pts col=upDist qcol="along+(upDist-length)" layer=2

Add new columns for the prediction points

bash v.db.addcolumn map=pred_pts layer=2 columns='pid INT, locID INT, ratio double precision, NEAR_X double precision, NEAR_Y double precision, elev double precision'

Now update the table

bash v.db.update pred_pts col=locID qcol="pid" layer=2 v.db.update pred_pts col=ratio qcol="round(along/length, 2)" layer=2 v.to.db map=pred_pts layer=2 option=coor col=NEAR_X,NEAR_Y v.what.rast map=pred_pts raster=lidar_dem column=elev layer=2

Use the v.db.univar function to find the maximum pid value.

bash v.db.univar -g map=snap_sites column=pid

Alternativitly, you can assign the values to the shell window to be used latter.

bash eval $(v.db.univar -g map=snap_sites column=pid)

This will give you all the stats from the map in the bash session and you can access them like:

bash echo $max

Now you can update the value of the pid columns which will be greater than all the values in the observations to make sure there are no duplicated pid values

bash v.db.update pred_pts col=pid qcol="$max+cat" layer=2

Step 7. Getting contributing area and covariate stats

to be continued…