Selecting Specific Elevation Areas Contiguous to the Coastline for Coastal Flood Mapping

Question

This is less a problem than an inquiry into a better way of doing something. I am curious how people might go about selecting, say, all areas less than or equal to 10 feet (from a DEM) that are contiguous to the coastline. The purpose of this exercise would be to map low-lying areas along the coast that may flood due to surge during a coastal storm.

Here is an abridged version of how I've been doing this sort of thing:

1) Extract from DEM all values less than or equal to 10 feet using the (raster) clip tool.
2) Convert the raster to a polygon using the Raster to Polygon tool.
3) Select areas contiguous to the coastline layer:
      a) Select all polygon areas that intersect the coast line. Export selection.
      b) Add to selection all areas that intersect the new selection. Export.
      c) Repeat this process until the difference in the number of attributes from one exported selection
          to the next reaches 0, meaning, you are no longer selecting new features when you select
          intersecting features.
4) Delete the remaining "orphaned" areas.

It's really step "3" that I'm trying to find an alternate method to. It can take up to 20 iterations of this step to select these areas for NYC. It just seems so tedious. I'm really curious how others might go about performing this task.

---

Update: I've come up with another method that eliminates the very iterative Step 3 from above. I am still curious to read how others might accomplish this task though.

Alternative Method:

1) Extract from DEM all values less than or equal to 10 feet using the (raster) clip tool.
2) Give every raster cell the same value (i.e. create a new field and have all values in that field = 1).
3) Convert raster to polygon using the Raster to Polygon tool.
4) Select and export features that intersect the coastline.

The above will get you your contiguous "flooded" area, but it will not have elevation values assigned to it, like the first step of steps above. To get those back, you'll need to do the following:

         A) Use the resulting features (shapefile or feature class) as a mask to clip the raster.
         B) Convert raster to polygon.

Answer 1

1. Use raster calculator Con(raster<30,1)

Apply Region Group in spatial analyst with this options:

• Convert output to polygons
• Select them by coastline
• Export selection to separate dataset and join it back to bigger set using Grid_Code field in both

Records that have match is what you need

Answer 2

This is going to require automation, and probably python is your best bet. Python has a while loop, that will repeat a set of commands as long as a condition is met. In this case, the condition would be no new selected features.

Here's an untested example. You'll need to assign paths to your feature classes. paste the paths between the r"...." at the start (Between Locals and End locals).

See comments (after hash (#)) for understanding of what's going on.

#--------------Locals

#Coastline feature class full path
coastLineFc = r"C:\test\test.gdb\Coastline"

#Polygon feature class from raster
elevFc = r"C:\test\test.gdb\DEMPolygon"

#Location for export
outLoc = r"C:\test\test.gdb"

#Output feature class name
#Feature class will have number extensions indicating each iteration
#i.e. output1, output2, etc etc
outName = "DEMPoly_iteration_"

#------------End locals

#Import your ArcGIS tools
import arcpy

#Make feature layer of elevation polygon feature class
#This allows you to select features
arcpy.MakeFeatureLayer_management (elevFc, "elevLyr")

#Initial selection - select polygons intersecting coast
arcpy.SelectLayerByLocation_management ("elevLyr", "", coastLineFc)

#Iteration variable - for naming feature class
i = 1

#Export first selection to new feature class
iOutName = outName + str(i)
arcpy.FeatureClassToFeatureClass_conversion ("elevLyr", outLoc, iOutName)
print "Created:", iOutName

#Counting varaibles
#Set initial count to zero
oldCount = 0
#Get new count from selection
newCount = int (arcpy.GetCount_management ("elevLyr").getOutput (0))

#While loop - select until no new selections
while oldCount != newCount:
    #Add 1 to iteration
    i += 1

    #Set old count to previous count
    oldCount = newCount

    #Select intersecting polygons from last selection
    arcpy.SelectLayerByLocation_management ("elevLyr", "", "elevLyr")

    #Check number of features selected
    newCount = int (arcpy.GetCount_management ("elevLyr").getOutput (0))

    #Check for new selection
    #Exit loop if no new selection
    #Makes while loop redundant but no harm/no foul
    if newCount == oldCount:
        break

    #Export (if new selection)
    iOutName = outName + str(i)
    arcpy.FeatureClassToFeatureClass_conversion ("elevLyr", outLoc, iOutName)
    print "Created:", iOutName

#Clean up
arcpy.Delete_management ("elevLyr")

#Done
print "Done"