Category Archives: GRASS

Creating Depth-Volume curves with GRASS-GIS

Our regional Drainage Authority prepared a reservoir at the mouth of a small dry riverbed to catch and regulate flood water coming from a mountain canyon. This reservoir was to act as a buffer to prevent flooding of agricultural fields and residential areas further down the valley. After a sudden rainstorm last week, the reservoir bravely fulfilled (pun intended ;-) ) it’s duty. Now we want to know how much water was actually captured, and to create a depth volume curve for the small “lake” that was formed. Here’s how I did this using GRASS.
Continue reading

Strahler stream order in GRASS

In hydrology, a stream network is composed of segments or “reaches” which are arranged in a hierachy. There are several systems of ordering the stream reaches, the most popular of which is the Strahler or Horton  number. GRASS GIS offers, alongside the watershed delineation tool r.watershed (discussed here), a set of addons for stream network analysis. We’ll examine how to use these addons, and how to use strahler ordering to improve the visual effect of a stream network map.

Continue reading

GIS on Scientific Linux 6

Exciting days these for open source GIS users. Debian 6 “Squeeze” is now the stable release, with full support for everything spatial. In the “rpm” camp, since Red Hat’s release of RHEL6 we’re looking forward to an equivalent selection of ready-to-go GIS software. The Fedora folks have been racing around the track to push out updated packages for each new Fedora release.
Developers of CentOS, the popular Red Hat clone, chose to give priority to security updates for the current CentOS 5.5 (a wise decision in my view). Although I’m sure that CentOS 6 is just a stone’s throw down the road. Meanwhile Scientific Linux has inched ahead with their version 6. Those looking to drive an enterprise class GIS workstation can choose the turn by turn steps below to setup GRASS GIS and Quantum GIS on SL 6.0.
So ladies and gents, rev up your engines and let’s get rolling.
Continue reading

Manipulating PostGIS layers from GRASS

GRASS users realize, right from the start, that GRASS maintains its own internal format for vector and raster layers. Whenever you want to use GRASS modules on GIS layers from some other source you first must import the data layers into GRASS, then run your analysis, and if necessary, export back to a required format.

But suppose you need to manipulate attributes of a vector layer that’s already stored in a PostGIS spatial database. A typical example: I have a GIS layer of rain gauges, stored in PostGIS. I need to add a column to the layer’s attribute table for elevation, and then get the elevation for each gauge from a GRASS based DEM (Digital Elevation Model) raster. Continue reading

Watershed Analysis with GRASS

Arcview users, needing to delineate watersheds and stream networks, would do well to choose the extension called  “Arc Hydro” (requires, at least the Spatial Analyst extension from ESRI). This extension introduces the concepts of “Batch Points” and “Adjoint Catchments”.  Batch points are locations that the user defines as drainage outlets. An adjoint catchment is the collection of all raster cells that drain into one of the batch points.

I’d like to demonstrate how to get similar results with GRASS GIS. As an example, we’ll create a set of catchments with their drainage outlets exactly at the points where the streams cross a road. We’ll assume our starting data includes an elevation raster called “dem” and a line vector called “roads”.  We first create the regular hydrology layers.

#set the region to the dem raster, and run the r.watershed module.
g.region -p rast=dem
r.watershed elev=dem drain=fdir basin=catch stream=str thresh=<your-threshold>

So far, pretty straightforward. There’s abundant information on r.watershed. I’ll just mention that the threshold value is the number of cells that will be the minimum catchment size. So if the resolution of our dem raster is, for example, 10×10 meters (each cell=100 sq. meters), then a threshold of 20,000 (=2,000,000 sq. meters) would create catchments of at least 2 sq. kilometers.

When the process finishes we’ll have three new raster maps: the flow direction map, the streams and the catchments. Let’s see what we’ve got so far:

#Convert the steams and catchments to vectors
r.to.vect in=catch out=catchments feature=area
# the stream raster usually requires thinning
r.thin in=str out=str_thin
r.to.vect in=str_thin out=streams feature=line
r.colors dem col=elevation
# Make a hillshade raster for displaying "3D"
r.shaded.relief map=dem shade=dem_shade zmult=1.5
# Now display layers
d.mon start=x0
d.his h=dem i=dem_shade
d.vect map=streams color=blue width=3
d.vect map=catchments type=boundary color=red
d.vect roads color=black width=2

Now we need to find all the points where streams cross roads. The v.overlay module does not deal with point vectors. Instead we use a trick in v.clean. When cleaning a line vector, all points where lines intersect and no node exists are considered “errors” and can be saved to a new point vector. So by merging the roads and streams vectors, we have lines (streams) crossing lines (roads) with no node at that location. Then we run v.clean to get all those intersection points in a new layer.

# Patch the streams and roads vectors together
v.patch in=streams,roads out=streams_roads
v.clean in=streams_roads out=streams_roads_clean tool=break error=cross_points
#View cross points on display
d.vect cross_points icon=basic/circle color=green size=12
# Save crossing points to a text file
v.out.ascii in=cross_points out=cross_points.txt format=point fs=space

Once we have the crossing points in a file, we simply run r.water.outlet in a loop to create the watersheds for each cross point.  The module r.water.outlet sets the raster value to 1 in the watershed and 0 elsewhere. We’ll use r.null to force the raster to ‘null’ outside of the watershed. THen we convert each raster to a vector (polygon).

i=0   #an iterator to give consecutive names to the new watersheds
while read X Y; do \
i=$(( ${i} + 1 ))  \
r.water.outlet drain=fdir east=$X north=$Y basin=wshed_$i  \
r.null -q wshed_$i setnull=0 \
r.to.vect -s in=wshed_$i out=wshed_$i feature=area \
g.remove rast=wshed_$i
done < cross_points.txt
echo "Created $i catchments"

We patch together all the area vectors and clean the merged watershed vector

v.patch in=`g.mlist vect pattern=wshed_* separator=,` out=wshed_patch
v.clean wshed_patch out=wshed_clean tool=snap,rmdupl thresh=20

Choose an appropriate snap threshold value based on your region resolution.  Additional manual cleaning may be required. After cleaning, reset the cat values (only for centroids):

v.category in=wshed_clean out=wshed_clean2 opt=del
v.category in=wshed_clean2 out=wshed_final type=centroid opt=add
g.remove vect=wshed_clean,wshed_clean2

Most likely we’ll want to calculate the area for each watershed.

# Add a table with a column for area in sq.km.
v.db.addtable wshed_final col="area_sqkm DOUBLE"
v.to.db map=wshed_final option=area col=area_sqkm unit=k

And finally, we can view the catchments, and their area values:

d.vect wshed_final type=boundary,centroid display=shape,attr attrcol=area_sqkm size=0 width=3 color=orange