How to add hypothetical values to a raster?

Question

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.

Answer 1

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:

1. Set environment extent to your DEM
2. Calculate raster called frood using DEM+OldDepth+1 (if it is all in feet)
3. Fill your DEM and calculate flow direction grid, call it FDIR
4. Modify grid paths in below script and run it.

Use raster calculator:

Con(“toriver”=>”DEM”,”toriver”-“DEM”)

To compute new depth:

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()