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
Pages: 1 2