I have 60 DEMs in a particular county created from aerial LiDAR that I'm trying to create watersheds for.One of the preliminary things I need to do is fill sinks, but first I need to know the fill limit for sinks in each DEM. To solve this problem, I read ESRI's section on creating depressionless DEMs and how to find sink depth. The problem is ESRI left me hanging as to which statistic is considerd the z limit.(http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z0000005q000000.htm).Below is the sinkdepth operation output for the first DEM. Which one of the following sinkdepth output statistics is the z limit value for Fill tool? <HistMin>1.1444091796875e-005</HistMin>
<HistMax>0.2707290649414063</HistMax>
<Approximate>0</Approximate>
<Metadata>
<MDI key="STATISTICS_MINIMUM">1.1444091796875e-005</MDI>
<MDI key="STATISTICS_MAXIMUM">0.27072906494141</MDI>
<MDI key="STATISTICS_MEAN">0.03058447859501</MDI>
<MDI key="STATISTICS_STDDEV">0.020851988436709</MDI>
The code below is what i used to create the sink depth raster including .xml statistics of each DEM. It takes all unfilled DEMs in a folder and outputs DEMs for sinkdepth:# ---------------------------------------------------------------------------
# find sink depth.py
# Created on: 2012-04-27 00:02:16.00000
# Author: T. Lupher
# Description:
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy, sys, os
from arcpy import env
from arcpy.sa import *
# Set workspace
env.workspace = "D:/DEM/LIDAR/"
# Allow overwrites
arcpy.env.overwriteOutput = True
# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")
# Create list of rasters
rasters = arcpy.ListRasters("*", "TIF")
rasterCount = len(rasters)
counter = 0
i = 1
try:
while counter < rasterCount:
for item in rasters:
# input/output directories
inputPath = "D:/DEM/"
flowdirFolder = "flowdirection/"
sinkFolder = "sinks/"
#pourPoints = "D:/DEM/pourpoints/alltraplocations.shp"
contribFolder = "sink_contributing_areas/"
zonestatFolder = "zonalstatistics/"
zonalfillFolder = "zonalfill/"
sinkdepthFolder = "sinkdepth/"
# input/output variables
outDropRaster = ""
flow = inputPath+flowdirFolder+"flowdir_"+str(i)+".tif"
sink = inputPath+sinkFolder+"sink_"+str(i)+".tif"
contrib = inputPath+contribFolder+"conarea_"+str(i)+".tif"
zstat = inputPath+zonestatFolder+"zonestat_"+str(i)+".tif"
zfill = inputPath+zonalfillFolder+"zonefill_"+str(i)+".tif"
sinkdepth = inputPath+sinkdepthFolder+"sinkdepth_"+str(i)+".tif"
# Process: Flow Direction
outFlow = FlowDirection(item, "NORMAL")
outFlow.save(flow)
# Process: Sink
outSink = Sink(flow)
outSink.save(sink)
# Process: Watershed
outWatershed = Watershed(flow, sink, "VALUE")
outWatershed.save(contrib)
# Process: Zonal Statistics
outZstat = ZonalStatistics(contrib, "VALUE", item, "MINIMUM", "DATA")
outZstat.save(zstat)
# Process: Zonal Fill
outZfill = ZonalFill(contrib, item)
outZstat.save(zfill)
# Process: Minus
outSinkDepth = Minus(outZfill, outZstat)
outSinkDepth.save(sinkdepth)
counter += 1
i += 1
except:
# If an error occurred while running a tool, then print the messages.
print arcpy.GetMessages()