25 June 2001
Keywords: CostDistance, flood, DEM, Avenue, Spatial Analyst, Map Calculator, map algebra
A reader writes,
I have a DEM of the Red River in Manitoba (5m cell, 15cm vert acc) and I would like to simulate a flood. I have a body of water fitting into the river and I would like it to rise in elevation and expand to fill all areas below and beside to simulate a growing body of water. When it breaches the river banks, I want it to expand and fill all areas with an elevation lower than it does, so long as it can reach it (eg not on the other side of a dam or dike).
There are many ways to do this in Spatial Analyst. One of the simpler uses just a few Map Calculator expressions, so you can test it out from the ArcView interface to see whether it works. I recommend trying it first on a small piece of the DEM (say, around 256 x 256 cells or so).
Begin with the DEM shown in a grid theme named "DEM". (Of course you can use other names, but I need a set of names to illustrate the process.) If the body of water is represented as a polygon theme, first convert that to a grid so that cells within the body have values and cells outside the body have NoData values. Let's call this theme "Source".
Suppose, for the sake of illustration, that the elevation of the body of water is 200 meters. Then the query
[DEM] > 200
executed in the Map Calculator produces a grid with the value 0 at elevations of 200 and lower and positive values (1) at higher elevations. Name the resulting theme "DEM > 200".
If [Source] and [DEM > 200] were polygon themes, we could select all contiguous points at elevations of 200 meters or less that are accessible to [Source] by using a simple theme-on-theme selection. Indeed, you could convert your themes to polygons and do this, but there is a better way using the grids directly.
The Map Calculator expression
([Source]).CostDistance([DEM > 200], NIL, NIL, NIL) = 0
will indicate the zone into which your source body will flood. If your source is defined consistently with the DEM, this will not augment the source at all, but it's probably a good place to start.
A brief explanation: this expression computes the cost of traversing the [DEM > 200] grid, beginning with any point in the [Source] grid. The cost is zero to travel through cells of elevation 200 meters or less, but otherwise the cost is 1. Therefore, the total travel cost will equal zero exactly within the region flooded by [Source].
Another way to describe this is that the CostDistance request is, fundamentally, a flooding request, but it accumulates a penalty for flooding across cells whose elevation exceeds 200 meters. Thus, by picking the cells of zero cost, you have found all cells reachable by crossing no elevations above 200 meters.
Iterate. Replace [Source] by the result of the previous map calculation. To accomplish this, you will need to reclassify the 1's as 0's and the 0's as NoData. The Log function does this, as in
Choose the next highest elevation for flooding: it depends on the map accuracy and precision. I will suppose, to illustrate, that a reasonable choice for the next highest level is 200.25 meters. Let's combine the two map calculations into one:
([Source]).CostDistance([DEM] > 200.25, NIL, NIL, NIL) = 0
(note the subtle difference in square brackets). The result shows the extent of flooding at the 200.25 meter elevation.
Keep iterating to the maximum flood elevation.
|220 meters--base level||229 meters--the floodway is inundated||233 meters--dikes separate the floodway from the spreading waters; islands appear in the flood|
(Click on any image to open a larger version in a new window.)
The process is readily automated with an Avenue script. I have appended the bare-bones version (no interface, little error checking, poor mechanism for user interruption, etc.).
A potential problem with this approach is that no flood maintains a constant elevation: it, too, is flowing and its elevation descends as the river descends. Therefore the calculation has to be considered of strictly local nature. You might be able to simulate a flood more realistically by adding a tiny positive value to the cost themes; then, you can limit the extent of the flood to n cells by selecting CostDistance values less than n times the tiny value. (Just make sure n * tiny is less than 1, so that the flood cannot hop over cells of higher elevation.) I have incorporated this in the script. If you don't want it, either remove the code or just set nTiny to zero.
There are more sophisticated methods one could use to make the model even more realistic, but I hesitate to go further: the best approach would be developed in conjunction with a water flow model.
CostDistance calculations can be time-consuming, so test cautiously: don't try to generate dozens of flood levels in one go. The process is iterative, so you can compute a flood to an intermediate elevation, inspect the results, and then resume from there.
Let me know if this helps.
The link below is a [Web page] showing an animation of the flood [and explaining how it was created]. The scenario shows what would happen if the Red River Floodway was unable to open (ice jam?) and divert floodwaters. http://www.geocities.com/cbouchar/interesting.html [about 4 MB].
That is a nice animation. It is very interesting to watch. Thank you for sharing it.
[Animation courtesy of Cameron Bouchard.]
Here is the Avenue script. At heart this script is a dozen lines--the kind of thing a GIS user routinely codes on the fly to get tasks done. It has been enhanced with comments, variables for modifying constant values, and a simple error check.
' Avenue script to generate successive flood levels.
' Requires Spatial Analyst.
' Modify the "model specification" values to suit the input.
' The source water body theme must be a grid theme named "Source".
' The elevation theme must be a grid theme named "DEM". Modify
' these names in this script as necessary.
' William A. Huber
' Quantitative Decisions
' 25 June 2001
theView = av.GetActiveDoc
' Obtain input grid themes.
thmSource = theView.FindTheme("Source") ' Should be NoData outside the
' ..source body
thmDEM = theView.FindTheme("DEM") ' Digital elevations.
' Obtain the model specification.
nSourceElevation = 200 ' The elevation corresponding to
' ..the source body
nFloodMax = 205 ' The maximum flood elevation;
' ..must not be less then nSourceElevation
nIncrement = 0.25 ' The flooding elevation increment
nExtent = 200 ' Do not flood outwards more than
' ..this many cells.
nTiny = 0.9/nExtent ' Used for limiting flood extents
' Perform the flooding.
gSource = thmSource.GetGrid
gDEM = thmDEM.GetGrid
for each xElevation in nSourceElevation..nFloodMax by nIncrement
' Flood to xElevation.
gFlood = (gSource.CostDistance((gDEM > xElevation)+nTiny, NIL, NIL, NIL)
if (gFlood.HasError) then
MsgBox.Info("Error at elevation" ++ xElevation.AsString,
' Display the flood stage.
thmFlood = gTheme.Make(gFlood)
thmFlood.SetName("Elevation" ++ xElevation.AsString)
' Prepare for the next step.
gSource = gFlood.log ' Converts 1 to 0, 0 to NoData
' end of script
26 June 2001: Created
23 July 2001: Updated link to flood animation