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.

Luckily, we had ordered in advance LIDAR coverage of the whole area, so I have a high resolution (1 meter per pixel) DTM of the reservoir and streams which feed into it during flash floods. I first visually examined the elevation values at the spillover point where the reservoir drains when it fills, and the lowest point in the center of the depression. I then ran the GRASS module r.lake feeding it the x,y coordinates of the deepest point as the “seed” and the elevation just below the spillover as the water level parameters. The result gave me the lake’s depths as a raster. Here’s the result:

amram reservoir

DTM and reservoir created with r.lake

Now I ran r.univar on the lake raster, to find

  • the maximum value = water depth
  • the sum of all pixels = total volume of the reservoir
  • the number of pixels = surface area covered by water

Since the region settings were 1 meter per pixel, the number of cells = area in sq.m. and the sum of all cells = volume in cubic meters.

Now I wanted to run the same procedure for water levels stepping down from the top spillover level to the bottom of the reservoir when it dries up due to leaching and evaporation. A script suggested itself, and here’s the depth_volume.sh script that I conjured up:


#!/bin/bash
# Author: micha, 7/11/2012; copyright: GPL >= 2
# Purpose: calculate volume, area and depth of a reservoir
# from a DEM raster and given water level using r.lake
# calculate volume and area with r.univar
# export the results to a text file for use in gnuplot
# Usage: depth_volume.sh {dem raster} {water level} {seed location as x,y}

if [ -z "$GISBASE" ] ; then
echo "You must be in GRASS GIS to run this program." 1>&2
exit 1
fi

if [ $# -lt 3 ] ; then
echo "Some inputs not defined. Usage: $0 dem_raster water_level seed_x,seed_y"
exit 1
fi

DEM=$1
LEVEL=$2
SEED=$3

g.message "Flooding lake at water level $LEVEL"
# Get resolution to calculate area and volume based on pixel size
NSRES=`g.region -g | grep nsres | cut -d= -f2`
EWRES=`g.region -g | grep ewres | cut -d= -f2`
PIXEL=`echo $NSRES*$EWRES | bc`
# Run r.lake
r.lake --quiet --o $DEM wl=$LEVEL xy=$SEED lake="$DEM"_tmp

# Collect needed numbers from r.univar
eval `r.univar -g "$DEM"_tmp`
VOLUME=`echo $sum*PIXEL/1000 | bc`
DEPTH=$max
AREA=`echo $n*$PIXEL | bc`

# Dump results into a text file
echo "$LEVEL $DEPTH $VOLUME $AREA" >> depth_volume.txt

The script takes three inputs: the DEM, water level and the X-Y coordinates for seeding the lake (as mentioned I chose the lowest spot in the depression). After running r.lake, the r.univar output is parsed to collect the needed numbers, and the results are dropped into a text file.
Now I ran this script in a loop with all the water levels that interested us. We wanted water volumes for each 1/2 meter drop in the water surface, so:

for level in 27 26.5 26 25.5 25 24.5 24 23.5 23 22.5; \
do depth_volume.sh dem $level 194365,397295; done

The depth_volume.txt file that is created has water level, volume and surface area for each of the above elevations. Now, to create a depth-volume curve, I could pull this text file into LibreOffice calc for designing a fancy chart. But for a quick and simple display, I used gnuplot with these commands:

set title "Depth-Volume Curve"
set xlabel "Volume 1000 m3"
set ylabel "Depth m."
set x2label "Area sqm."
set key left box
set x2tics nomirror
set autoscale x2
set grid y
set terminal png enhanced size 800,600
set output "depth_volume.png"
plot "depth_volume.txt" using 3:2 title "Volume 1000 m3" with lines lw 3, \
"depth_volume.txt" using 4:2 title "Area m2" with lines axes x2y1

The graph looks like this:

depth_volume

Depth Volume curve from gnuplot

Once I have setup the GRASS script and gnuplot command file, running the model with other parameters or in some other location is a piece of cake.

Creative Commons License
This work, unless otherwise expressly stated, is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
  1. #1 by Nimrod Cnaan on 28/11/2012 - 22:55

    Very impressive

  2. #2 by Micha Silver on 28/11/2012 - 22:57

    Thanks!

  3. #3 by dbb on 29/11/2012 - 00:08

    great recipe Micha !

  4. #4 by Marcello Benigno on 20/12/2012 - 01:29

    Micha very good, thanks once again for sharing your knowledge!

(will not be published)