I have to create Non-Regulatory Products for a few FEMA regions. Right now I am working on the 'Increased Flooding Scenarios' where you add a hypothetical value to the existing Base Flood (1', 2', and 3', specifically). I already created the Base Flood (or Coastal Depth Grid), now I am stuck on how to add a hypothetical value to it (1 foot) and have a 'plus 1' raster extending to the right elevation. I've tried various methods where the elevation increases by a foot, but the raster does not extends to the right elevation (Con, Plus, Raster Calculator). Attached is a flowchart in FEMA's guidelines on how to create the flooding scenarios.
-
Please tag with your software and version. – Tom Mar 2 '16 at 18:02
I assume that inundation elevation is different along the coast:
Otherwise it is simply finding areas where elevation below new level.
There are multiple solutions, e.g. calculate focal statistics (mean), mosaic it with existing flood level (top) and repeat until no new DEM cells are below mosaic.
I use my own:
- Set environment extent to your DEM
- Calculate raster called frood using DEM+OldDepth+1 (if it is all in feet)
- Fill your DEM and calculate flow direction grid, call it FDIR
- Modify grid paths in below script and run it.
Use raster calculator:
Con(“toriver”=>”DEM”,”toriver”-“DEM”)
import arcpy, os, traceback, sys, numpy
from arcpy import env
from arcpy.sa import *
fdir=r'D:\Scratch\fdir'
demG=r'D:\Scratch\flood'
d=arcpy.Describe(demG)
try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
dirArray = arcpy.RasterToNumPyArray(fdir,"","","",-9999)
dem = arcpy.RasterToNumPyArray(demG,"","","",-9999)
nRows,nCols=dirArray.shape
blankArray=arcpy.RasterToNumPyArray(demG,"","","",-9999)
cellsTotal=nCols*nRows
d=arcpy.Describe(fdir)
origin=d.extent.lowerLeft
cSize=arcpy.Raster(fdir).meanCellHeight
# directions to find cell neighbour
fDirs=(1,2,4,8,16,32,64,128)
dCol=(1, 1, 0, -1, -1,-1, 0,1)
dRow=(0, -1, -1, -1, 0, 1, 1,1)
# flipped
dRow=(0, 1, 1, 1, 0, -1, -1,-1)
aDict={}
for i,v in enumerate(fDirs):
aDict[v]=(dCol[i],dRow[i])
# main loop
arcpy.SetProgressor("step", "", 0, nRows)
for nRow in range (nRows):
for nCol in range (nCols):
S=dirArray[nRow,nCol]
if S in (-1,-9999): continue
nR,nC=nRow,nCol
cells=[(nR,nC)]
z=dem[nR,nC]
while True:
direction=dirArray[nR,nC]
dX,dY=aDict[direction];nC+=dX
if nC not in range(nCols): break
nR+=dY
if nR not in range(nRows): break
cells.append((nR,nC))
S=dirArray[nR,nC]
z=dem[nR,nC]
if S<-9998 or z>-9998:
break
if S==-1:
z=blankArray[nR,nC]
break
for nR,nC in cells:
blankArray[nR,nC]=z
dirArray[nR,nC]=-1
arcpy.SetProgressorPosition()
myRaster = arcpy.NumPyArrayToRaster(blankArray,origin,cSize,cSize)
oneGrid=Con(myRaster<>-9999,myRaster)
oneGrid.save(r'D:\Aerials\nztm\toriver')
del dirArray,blankArray
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
-
thanks for your response! Digging on the internet I found a methodology that fits exactly what the client wants. floods.org/Files/Conf2014_ppts/D6_Donnelly.pdf – Naomy P. Mar 8 '16 at 18:08