Predict Discharge With Unit Hydrographs Automated

import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

def Isochrone(Slope, FlowDirection, FlowAccumulation, Point, SnapDist, MaxTime, TimeInterval, Inter_Data):

    # Returns isochrone raster
    # Slope = Slope raster (% rise)
    # FlowDirection = Flow direction raster
    # FlowAccumulation = Flow accumulation raster
    # Point = Pour point
    # SnapDist = snap distance
    # Inter_Data = 'False' = Delete intermediate data

    def PrintLastMsg(tool):

        # Prints a tool string passed to it then arcpy messages
        # tool = string

        count = arcpy.GetMessageCount()
        print('')
        print(tool)
        print(arcpy.GetMessage(1))
        print(arcpy.GetMessage(count-1))


    # Snap the pour point to the highest flow accumulation within snap distance

    pointSnap = arcpy.sa.SnapPourPoint(Point, FlowAccumulation, SnapDist)
    PrintLastMsg('SnapPourPoint')

    # Delineate the pour point watershed

    pt_Watershed = arcpy.sa.Watershed(FlowDirection, pointSnap, 'VALUE')
    PrintLastMsg('Watershed')
    pt_Watershed.save('pt_Watershed')
    pointSnap.save('pointSnap')

    # Apply a mask to make raster processing faster for intermediate data
    arcpy.env.mask = pt_Watershed

    # -----------------------------------------------------------------
    # Perform slope area term calculation for culvert specific watershed

    SlopeAreaTerm = (SquareRoot(Slope) * SquareRoot(FlowAccumulation))
    SlopeAreaTerm.save('SlopeAreaTerm')

    # Find the average slope of watershed

    SlopeAnalysis = arcpy.GetRasterProperties_management(SlopeAreaTerm, "MEAN")
    PrintLastMsg('GetRasterProperties_management')
    SlopeMean = SlopeAnalysis.getOutput(0)
    SlopeMeanFt = float(SlopeMean)
    print(SlopeMeanFt)

    # Create unlimited velocity feature (Results in
    # unrealistic values so further analysis is required

    VelocityUnlimited = 0.1 * (SlopeAreaTerm / SlopeMeanFt)
    VelocityUnlimited.save('VelocityUnlimited')

    # Limit the maximum velocity

    VelocityLim = Con(VelocityUnlimited, VelocityUnlimited, "0.02", '[Value] >= 0.02')
    PrintLastMsg('Con')
    VelocityLim.save('VelocityLim')

    Velocity = Con(VelocityLim, VelocityLim, "2", '[Value] < 2')
    PrintLastMsg('Con')
    Velocity.save('Velocity')

    # Create Flow Length raster. This determines how long each 
    # cell in the watershed takes to reach the pour point


    WeightRaster = 1 / Velocity
    FlowLen = arcpy.sa.FlowLength(FlowDirection, "DOWNSTREAM", WeightRaster)
    PrintLastMsg('FlowLength')
    WeightRaster.save('WeightRaster')
    FlowLen.save('FlowLengthRaster')

    # -----------------------------------------------------------------
    # Create a remap table with intervals of TimeInterval until MaxTime is reached

    RemapList = []
    countList = 0

    while countList < MaxTime/TimeInterval:
        RemapList.append(
            [
                float(countList*TimeInterval), 
                float((countList+1)*TimeInterval), 
                (countList+1)*TimeInterval
            ]
        )
        countList = countList + 1

    timeRemap = RemapRange(RemapList)

    # Reclassify Flow Length raster to organize into classes.

    Isochrone = arcpy.sa.Reclassify(FlowLen, 'VALUE', timeRemap, 'NODATA')
    PrintLastMsg('Reclassify')

    # Add fields, calculate flow area per 10 minute increment 
    # and average flow per second (m^3) for each class

    desc = arcpy.Describe(Isochrone)
    CellArea = desc.meanCellWidth * desc.meanCellWidth

    arcpy.AddField_management(Isochrone, "HydroOrdin", "DOUBLE")
    arcpy.CalculateField_management(
        Isochrone, 
        "HydroOrdin", 
        "!COUNT! * {0} / {1}".format(CellArea, TimeInterval), 
        "PYTHON")

    # Delete intermediate data

    if Inter_Data == 'False':

        print('Deleting Intermediate Data')
        arcpy.Delete_management('pointSnap')
        arcpy.Delete_management('pt_Watershed')
        arcpy.Delete_management('SlopeAreaTerm')
        arcpy.Delete_management('VelocityUnlimited')
        arcpy.Delete_management('VelocityLim')
        arcpy.Delete_management('Velocity')
        arcpy.Delete_management('WeightRaster')
        arcpy.Delete_management('FlowLengthRaster')

    return Isochrone

Leave a comment