r.damflood.md

r.damflood example

I am trying to run r.damflood. Unfortunately the documentation for this addon has no examples and I can't figure out what I need to use as input.

According to the documentation you must prepare your input using the r.dam addon. The problem is that r.dam only exists on GRASS 6, while r.damflood only exists on GRASS 7. Fortunately, r.dam code is not really complicated so you can convert it to GRASS 7 compatible code with relateive ease.

This is what I've tried to do in this notebook.

Note: I have changed most of the map names in order to translate them to english from italian.

r.dam replication

Set up

1. Open the Northa Carolina LOCATION
2. Set the region:

g.region -p raster=elevation

1. Create the map containing the manning coefficients. For simplicity, we are going to use a constant value throughout the map:

r.mapcalc --overwrite exp="$manning=if(elevation > 0, 0.1)" 

Input

Create the following variables:

units=meters            # input
resolution=10.          # input

damlevel=105            # input
waterlevel=100          # input
dtm=elevation           # input
manning=manning         # input

# these are the coordinates of a point within the reservoir
# If you choose to do the dam in a different place than the
# one indicated in this example you will have to change them
EC=637798.742138,       # input
NC=218448.113208        # input

total_dam=total_dam     # output

dam_v=dam_v             # input
dam_r=dam_r
dam_buffer=dam_buffer
dam_dtm=dam_dtm

dambreak_v=dambreak_v
dambreak_r=dambreak_r
dambreak_dtm=dambreak_dtm

lake_dtm=lake_dtm
lake_map=lake_map
lake_map_no_buffer=lake_map_no_buffer
lake_volume=lake_volume

Procedure

1. Go to the digitizer and create the dam_v vector which describes the dam geometry

2. Convert the vector representation of the dam to a raster. We also create a DTM which contains the dam geometry and a buffer (I am not sure why we need the buffer but this is what r.dam was doing):

v.to.rast --overwrite --quiet input=$dam_v output=$dam_r use=cat type=point,line
r.mapcalc --overwrite exp="$dam_r=if(isnull($dam_r), 0, $dam_r)"
r.mapcalc --overwrite exp="$dam_dtm=if($dam_r>0 && $dtm<$damlevel, $damlevel, $dtm)"
r.buffer -z --overwrite --quiet input=$dam_r output=$dam_buffer distances=$resolution units=$units
r.mapcalc --overwrite exp="$dam_buffer=if(isnull($dam_buffer), 0, $dam_buffer)"

The dam_dtm will be something like this:

While the dam_buffer will be something like this:

1. Create the lake_dtm.

r.mapcalc --overwrite exp="$lake_dtm=if($dam_buffer>0 && $dtm<$damlevel, $damlevel, $dtm)"

AFAIK the only difference between this one and the dam_dtm is the usage of the buffered dam:

1. Calculate the lake_map:

r.lake elevation=$lake_dtm wl=$waterlevel lake=$lake_map coordinates=$EC,$NC --overwrite
r.mapcalc --overwrite exp="$lake_map=if(isnull($lake_map), 0, $lake_map)"

r.lake gives the following output:

Lake depth from 0.000000 to 16.067245 (specified water level is taken as zero)
Lake area 1433000.000000 square meters
Lake volume 7868751.045227 cubic meters

While lake_map looks like this (we also display the dam for reference):

The following steps are taken directly from r.dam. (AFAIK lake_map_no_buffer is being used in order to calculate the "correct" volume while the total_dam is the height of the dam below the water. I don't think that these are useful for r.damflood).

r.mapcalc --overwrite exp="$lake_map_no_buffer=if($dam_dtm<$waterlevel && ($lake_map[-1,0]>0 || $lake_map[0,1]>0 || $lake_map[1,0]>0 || $lake_map[0,-1]>0), $waterlevel - $dam_dtm, $lake_map)"
r.mapcalc --overwrite exp="$lake_volume=$lake_map_no_buffer * ewres() * nsres()"
r.mapcalc --overwrite exp="$total_dam=$dam_dtm - $dtm"
r.univar $lake_volume   # the "sum" attribute is the total volume of the lake.

The "correct" lake volume is: 7900275.46234131

total_dam:

At this point, r.dam goes on to delete the following maps:

• dam_buffer
• lake_map
• lake_volume
• dam_r
• lake_dtm

I.e from this point on we can only use these maps:

• dam_dtm
• lake_map_no_buffer
• total_dam

Dambreak

r.damflood needs a map defining the "dambreak". I am not really sure how I should create this. I checked the source code and it seems to me that "dambreak" should be an elevation like map but I am not sure about that. Anyway, this is what I did:

1. Create a copy of dam_v using:

v.split --overwrite input=dam_v output=dambreak_v length=50 units=meters

1. Go to the digitizer and edit dambreak_v. Remove the segments where the dam will get breached:

1. Convert the vector representation of the dambreak to a raster + a DTM (this step was NOT being done in r.dam):

v.to.rast --overwrite input=$dambreak_v output=$dambreak_r use=cat type=point,line
r.mapcalc --overwrite exp="$dambreak_r=if(isnull($dambreak_r), 0, $dambreak_r)"
r.mapcalc --overwrite exp="$dambreak_dtm=if($dambreak_r>0 && $dtm<$damlevel, $damlevel, $dtm)"

The dambreak_dtm will be like this:

r.damflood

We can now go on to run r.damflood. I've tried several things, this seems to run but the results are obviously wrong... (the water is not flowing from the dam breach!):

r.damflood elev=dam_dtm lake=lake_map_no_buffer manning=manning dambreak=total_dam tstop=40 deltat=1 h=wh_ vel=wv_ hmax=whmax vmax=wvmax imax=wimax wavefront=wavefront method='dambreak-without_hypotesis' --quiet --overwrite

1. t = 1s:
2. t = 10s:
3. t = 20s:
4. t = 40s: