Community

Help with Script to Loop through DEMs and get Raster Properties

2908
1
10-04-2012 05:34 AM
KateLyndegaard
by
New Contributor III
‎10-04-2012 05:34 AM 
#01/10/2012
#Script to automate terrain analysis
#Kate Lyndegaard

#++++++++++++++++++++++++++++++++++++++++++++++++

import arcpy
import os
import traceback
import sys
import fnmatch

from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

#DEM
demPath = arcpy.GetParameterAsText(0)

#Output folder
outputFolder = arcpy.GetParameterAsText(2)

#Set workspace
arcpy.env.workspace = arcpy.GetParameterAsText(1)

try:
    # Execute Int
    outDEM = Int(demPath)

    # Save the output 
    outDEM.save("outDEM")    

    #Split raster
    arcpy.SplitRaster_management(outDEM, outputFolder, "Tile", "SIZE_OF_TILE", "GRID", "BILINEAR", "", "5280 5280","","FEET")    
    
except:
    # Get the traceback object
    #
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]

    # Concatenate information together concerning the error into a message string
    #
    pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
    msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"

    # Return python error messages for use in script tool or Python Window
    #
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)

try:
    #Output folder
    outputFolder = arcpy.GetParameterAsText(2)

    #Set workspace
    arcpy.env.workspace = outputFolder

    demList = arcpy.ListRasters()

    for dem in demList:
        elevSTDresults = arcpy.GetRasterProperties_management(dem,"STD")
        elevSTD = elevSTDresults.getOutput(0)
        elevSTDNUM = float(elevSTD)

        if elevSTDNUM < 66:
            arcpy.env.workspace = outputFolder
            arcpy.Rename_management(dem, "flat_" + dem)
        elif elevSTDNUM > 66 and elevSTDNUM < 328:
            arcpy.env.workspace = outputFolder
            arcpy.Rename_management(dem, "phill_" + dem)
        elif elevSTDNUM > 328:
            arcpy.env.workspace = outputFolder
            arcpy.Rename_management(dem, "hill_" + dem)except:
    
    # Get the traceback object
    #
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]

    # Concatenate information together concerning the error into a message string
    #
    pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
    msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"

    # Return python error messages for use in script tool or Python Window
    #
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)           
       
try:
    #Output folder
    outputFolder = arcpy.GetParameterAsText(2)

    #Set workspace
    arcpy.env.workspace = outputFolder
    
    demList1 = arcpy.ListRasters()
    
    for dem in demList1:
        if fnmatch.fnmatch(dem,"flat_*"):
            arcpy.env.workspace = outputFolder
            arcpy.RasterToPolygon_conversion(outputFolder + "\\" + dem, outputFolder+ "\\" + dem +".shp", "SIMPLIFY")
        elif fnmatch.fnmatch(dem,"phill_*"):
            arcpy.env.workspace = outputFolder
            arcpy.RasterToPolygon_conversion(outputFolder + "\\" + dem, outputFolder+ "\\" + dem +".shp", "SIMPLIFY")
        elif fnmatch.fnmatch(dem,"hill_*"):
            arcpy.env.workspace = outputFolder
            arcpy.RasterToPolygon_conversion(outputFolder + "\\" + dem, outputFolder+ "\\" + dem +".shp", "SIMPLIFY")

except:
    # Get the traceback object
    #
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]

    # Concatenate information together concerning the error into a message string
    #
    pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
    msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"

    # Return python error messages for use in script tool or Python Window
    #
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)

try:
    #Output folder
    outputFolder = arcpy.GetParameterAsText(2)
    
    #Set workspace
    arcpy.env.workspace = outputFolder

    #Create folders to hold merged output     
    flatFolder = arcpy.CreateFolder_management(outputFolder, "Flat")    
    partFolder = arcpy.CreateFolder_management(outputFolder, "Partly_Hilly")
    hillFolder = arcpy.CreateFolder_management(outputFolder, "Hilly")
    
    flatFile = arcpy.ListFeatureClasses("flat*", "Polygon")
    tempListf = []
    for list in flatFile:
        item = list + ';'
        tempListf.append(item)
        finalListf = ''.join(tempListf)
        mergeStrF = finalListf.rstrip(';')
        mergeStrF2 = outputFolder + "\\" + mergeStrF
    
    phillFile = arcpy.ListFeatureClasses("phill*", "Polygon")
    tempListp = []
    for list in phillFile:
        item = list + ';'
        tempListp.append(item)
        finalListp = ''.join(tempListp)
        mergeStrP = finalListp.rstrip(';')
        mergeStrP2 = outputFolder + "\\" + mergeStrP
        
    hillFile = arcpy.ListFeatureClasses("hill*", "Polygon")
    tempListh = []
    for list in hillFile:
        item = list + ';'
        tempListh.append(item)
        finalListh = ''.join(tempListh)
        mergeStrH = finalListh.rstrip(';')
        mergeStrH2 = outputFolder + "\\" + mergeStrH
        
    # Provide name and path of the merged feature class    
    flatMerged = outputFolder + "\\" + "Flat" + "\\" + "flatMerged.shp"
    phillMerged = outputFolder + "\\" + "Partly_Hilly" + "\\" + "phillMerged.shp"
    hillMerged =outputFolder + "\\" + "Hilly" + "\\" + "hillMerged.shp"

    if not tempListf:
        pass
    else:
        arcpy.Merge_management(mergeStrF2, flatMerged)

    if not tempListp:
        pass
    else:
        arcpy.Merge_management(mergeStrP2, phillMerged)

    if not tempListh:
        pass
    else:
        arcpy.Merge_management(mergeStrH2, hillMerged)
              
except:

    # Get the traceback object
    #
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]

    # Concatenate information together concerning the error into a message string
    #
    pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
    msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"

    # Return python error messages for use in script tool or Python Window
    #
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)


After running the above code, I'm getting the error: ExecuteError: ERROR 001100: Failed because no statistics is available.
Failed to execute (GetRasterProperties).

I'm not sure why I'm getting this error as statistics have been calculated for the input DEM. The script seems only to read the first few tiles of the split raster output- then it fails. Does anyone have any suggestions?

Thank you
Tags (2)
0 Kudos
1 Reply
VandanaRaghunathan
by
New Contributor III
‎03-12-2015 09:00 AM

Hello Kate,

Have you tried using Calculate Statistics tool before running the GetRasterProperties?

for dem in demlist:
    arcpy.CalculateStatistics_management(dem)
     elevSTDresults = arcpy.GetRasterProperties_management(dem,"STD")

Thanks,

Vandana

0 Kudos