Depressionless DEMs: which statistic is considered the z Limit value for Fill tool?

274
0
04-28-2012 02:08 PM
BL
by
New Contributor II
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()  
0 Kudos
0 Replies