GISHydroNXT System Documentation
#*********************************************************************************
# Author: UMD
# Date: 24-07-2018
# Modified: n/a
# Classes: CalculateAttributes()
# Functions: hydro.expandExtent() ; hydro.TcFC() ; hydro.TcST_MRLC()
# hydro.TcST() ; hydro.TcImpMRLC() ; hydro.TcImpUltimate()
# hydro.TcImpUSGS() ; hydro.TcImpAnderson()
# Modules: arcpy ; os ; time
# Comments: n/a
#*********************************************************************************
class CalculateAttributes(object):
"""Implementation for GISHydroNXT_addin.button10 (Button)"""
def __init__(self):
self.enabled = False
self.checked = False
def onClick(self):
# For some reason the workspace must be setted again
arcpy.env.scratchWorkspace = scratchfolder
arcpy.env.workspace = optfolder
#*******************************************************************************************************
# input data needed for "subrivers.shp" and "subsheds.shp"
#*******************************************************************************************************
dem = optfolder + "/dem"
subrivers = optfolder + "/subrivers.shp"
slope_stats = optfolder + "/slope_stats.dbf"
polyras = optfolder + "/polyras"
elevzones = optfolder + "/elevzones.shp"
elevmerge = optfolder + "/elevmerge.shp"
flowdir = optfolder + "/flowdir"
fc = optfolder + "/subshed.shp" # shapefile
subsheds = optfolder + "/subsheds" # raster data
outlets = optfolder + "/outlets"
longfp = optfolder + "/longfp.dbf"
slope_sheds = optfolder + "/slope_sheds.dbf"
slopegrid = optfolder + "/slope_calc/landslope"
cntable = optfolder + "/cntable.dbf"
lu = optfolder + "/landuse"
#*******************************************************************************************************
# add length, and slope fields. Update "slope" field
# "subrivers.shp" attributes calculation
#*******************************************************************************************************
arcpy.AddField_management(subrivers,"Length","FLOAT",15,4)
arcpy.AddField_management(subrivers,"Slope","Double",15,9) # 08/05/2013: "Float" to "Double" and field scale from 4 to 9
arcpy.CalculateField_management(subrivers,"Length","!shape.length@feet!","PYTHON") # length converted into feet
arcpy.PolylineToRaster_conversion(subrivers,"ARCID",polyras,"MAXIMUM_LENGTH","NONE",30)
arcpy.RasterToPolygon_conversion(polyras,elevzones,"NO_SIMPLIFY","Value")
arcpy.Dissolve_management(elevzones, elevmerge, "GRIDCODE") # 08/05/2013: added to handle multiple polygons
arcpy.sa.ZonalStatisticsAsTable(elevmerge,"FID",dem,slope_stats,"DATA","ALL") # NODATA is changed to DATA # 08/05/2013: elevzones changed to elevmerge
zonalcur = arcpy.SearchCursor(slope_stats,"","","MIN;MAX","")
lencursor = arcpy.arcpy.SearchCursor(subrivers,"","","Length","")
minlst = []
maxlst = []
global lenlst
lenlst = []
for z in zonalcur:
min_elev = (z.getValue("MIN")) # removed multiplication factor of 3.2808 because
max_elev = (z.getValue("MAX")) # elevation units are already in feets
minlst.append(min_elev)
maxlst.append(max_elev)
for l in lencursor:
len_cur = l.getValue("Length")
lenlst.append(len_cur)
diff = [maxlst - minlst for maxlst, minlst in zip(maxlst, minlst)]
global slopeval
slopeval = [diff/lenlst for diff, lenlst in zip(diff, lenlst)]
global no_subwatersheds
no_subwatersheds = len(slopeval)
# slope values should at least be 0.0001
for index, item in enumerate(slopeval):
if item == "0":
slopeval[index] = "0.0001"
count = 0
fc = optfolder + "/subrivers.shp"
slopecur = arcpy.UpdateCursor(fc)
for s in slopecur:
val = slopeval[count]
s.Slope = val
slopecur.updateRow(s)
count = count + 1
# update "GRID_CODE" using "FROM_NODE" attribute values
nodes = arcpy.UpdateCursor(subrivers)
for node in nodes:
node.setValue("GRID_CODE", node.getValue("FROM_NODE"))
nodes.updateRow(node)
del node
del nodes
#*******************************************************************************************************
# Fields added:
# 1] length
# 2] slope
# 3] Tc method
# 4] Area in square miles
# 5] CN
#
# Update "slope" field "subsheds.shp" attributes calculation
#*******************************************************************************************************
flds = arcpy.sa.FlowLength(flowdir, "DOWNSTREAM","#")
minflds = arcpy.sa.ZonalStatistics(subsheds, "VALUE", flds, "MINIMUM", "DATA")
fldswo = arcpy.sa.Minus(flds, minflds)
NotOutlet = arcpy.sa.IsNull(outlets)
fdrnooutlet = arcpy.sa.Divide(flowdir, NotOutlet)
fluswb = arcpy.sa.FlowLength(fdrnooutlet, "UPSTREAM","#")
flplus = arcpy.sa.Plus(fldswo, fluswb)
lengthlfp = arcpy.sa.ZonalStatistics(subsheds, "VALUE", flplus, "MAXIMUM", "DATA")
length_lfp = lengthlfp* 3.2808
# zonal statistics for "longfp" and "slope_sheds"
slopegrid = arcpy.sa.Divide(slopegrid, 3.2808) # 07-06-2015: Landslope value should be factored into meters unit
arcpy.sa.ZonalStatisticsAsTable(subsheds, "VALUE", length_lfp, longfp, "DATA", "ALL") # 07-25-2013: "NODATA" was changed to "DATA"
arcpy.sa.ZonalStatisticsAsTable(subsheds, "VALUE", slopegrid, slope_sheds, "DATA", "ALL") # 07-25-2013: "NODATA" was changed to "DATA"
# calculate area in square miles
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"AreaMi2","FLOAT",15,4)
arcpy.CalculateField_management(subshed_poly,"AreaMi2","!shape.area@squaremiles!","PYTHON")
# update TcMethod
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"TcMethod","TEXT",15,4)
subpoly = arcpy.UpdateCursor(subshed_poly)
tc_name_lst = []
tc_name_lst.append(Tc_method)
tc_records = len(slopeval) # 09-30-2013: Use slopeval list to get length of records to insert in Tc Method list
tc_n_list = tc_name_lst * tc_records
tc_index = 0
for tc in subpoly:
Method = tc_n_list[tc_index]
tc.TcMethod = Method
subpoly.updateRow(tc)
tc_index = tc_index + 1
# update CN
cngrid = optfolder + "/curvenumber"
arcpy.sa.ZonalStatisticsAsTable(subsheds, "VALUE", cngrid, cntable, "NODATA", "ALL")
cntable = optfolder + "/cntable.dbf"
cn_mean = arcpy.SearchCursor(cntable,"","","MEAN","")
global cn_list
cn_list = []
for cn in cn_mean:
mean_cn = cn.getValue("MEAN")
cn_list.append(mean_cn)
cn_lst = 0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"CurveNum","FLOAT",15,4)
cn_update = arcpy.UpdateCursor(subshed_poly)
for cnum in cn_update:
sheds_cn = cn_list[cn_lst]
cnum.CurveNum = sheds_cn
cn_update.updateRow(cnum)
cn_lst = cn_lst + 1
# update slope -- get mean values from "slope_sheds", add field in "subshed_poly", and update it
slope_sheds = optfolder + "/slope_sheds.dbf"
slope_mean = arcpy.SearchCursor(slope_sheds,"","","MEAN","")
slope_list = []
for slp in slope_mean:
mean_slope = slp.getValue("MEAN")
slope_list.append(mean_slope)
# assume a minimum slope condition [MDSHA.TimeOfConcentration.CalculateTc -- line 44]
slope_list = [0.001 if i == 0 else i for i in slope_list]
# adjust units to "feet/feet"
slope_list = [i/3.2808 for i in slope_list]
# create a list holding subshed (raster) count fields
subarea = []
rows = arcpy.SearchCursor(subsheds,"","","Count","")
for row in rows:
thecount = row.getValue("Count")
subarea.append(thecount)
#*******************************************************************************************************
# calculate Tc at the end as we need Slope, CN, and LngFlwPth (L) for its computation
#*******************************************************************************************************
# SCS method
if Tc_method == "SCS Lag Formula":
# update LFP
longfp = optfolder + "/longfp.dbf"
lfp_max = arcpy.SearchCursor(longfp,"","","MAX","") # already in feets
lfp_list = []
for lfp in lfp_max:
max_lfp = lfp.getValue("MAX")
lfp_list.append(max_lfp)
lfp_lst = 0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"LngFlwPth","FLOAT",15,4)
lfp_update = arcpy.UpdateCursor(subshed_poly)
for lonfp in lfp_update:
sheds_lfp = lfp_list[lfp_lst]
lonfp.LngFlwPth = sheds_lfp
lfp_update.updateRow(lonfp)
lfp_lst = lfp_lst + 1
slp_lst = 0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"Slope","FLOAT","#","#")
slope_update = arcpy.UpdateCursor(subshed_poly)
for slope in slope_update:
sheds_slope = slope_list[slp_lst]
slope.Slope = sheds_slope
slope_update.updateRow(slope)
slp_lst = slp_lst + 1
global tc_list
tc_list = []
for a,b,c in zip(lfp_list, cn_list, slope_list):
tc_val = (np.float64(100*((a)**0.8)*(((1000/b)-9)**0.7))/(19000*((c)**0.5)))/60 # convert float to np.float to handle zero division error
tc_list.append(tc_val)
tc_lst=0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"Tc","FLOAT",15,4)
tc_update = arcpy.UpdateCursor(subshed_poly)
for tc in tc_update:
sheds_tc = tc_list[tc_lst]
tc.Tc = sheds_tc
tc_update.updateRow(tc)
tc_lst = tc_lst + 1
# Hydrology Panel Method
elif Tc_method == "Hydrology Panel Tc Method":
# get "gc_lst" before extent definition to avoid intersection for only sub-basin
Mdprov = r"" + Directory + "/data/maryland/Mdprov.shp"
subshed = optfolder + "/subshed.shp"
subshed_prov = optfolder + "/subshed_prov.shp"
arcpy.Intersect_analysis([Mdprov,subshed],subshed_prov,"ALL","#","INPUT")
arcpy.AddField_management(subshed_prov,"Area","FLOAT",15,4)
arcpy.CalculateField_management(subshed_prov,"Area","!shape.area@squaremeters!","PYTHON")
# clip and save sub-sheds 10% bigger than the original extent
fc = optfolder + "/subshed.shp"
desc = arcpy.Describe(fc)
rows = arcpy.SearchCursor(fc)
FC = []
ST = []
IA = []
lst = 0
for idx,row in enumerate(rows):
aPoly = row.getValue(desc.shapefieldname)
setExtent = aPoly.extent
expPercent = 0.1
bigExtent = hydro.expandExtent(setExtent,expPercent)
arcpy.env.extent = arcpy.Extent(int(bigExtent[0]), int(bigExtent[1]), int(bigExtent[2]), int(bigExtent[3]))
facc = optfolder + "/flowacc"
faccdesc = arcpy.Describe(facc)
cellSize = faccdesc.meanCellHeight
arcpy.env.cellSize = cellSize
# create mask raster for each sub-basin [04-18-2014: Legacy does same masking and only allow calculations for sub-basin extent]
basingrid = optfolder + "/basingrid"
arcpy.MakeFeatureLayer_management(fc, "layer" + str(idx), ' "FID" = ' + str(idx))
mask = arcpy.Clip_management(basingrid, "#", optfolder + "/shedMask" + str(idx), "layer" + str(idx), "0", "ClippingGeometry")
# clip and save sub-sheds 10% bigger than the original extent
arcpy.Clip_management(lu,"%f %f %f %f" %(bigExtent[0], bigExtent[1], bigExtent[2], bigExtent[3]), optfolder + "/Tc_temp"+str(idx), "#", "#", "NONE")
Tc_sub = arcpy.sa.Times(optfolder + "/Tc_temp"+str(idx),mask)
Tc_sub.save(optfolder + "/Tc_subshed"+str(idx))
# Determine FC, ST, Impervious counts, and save them to separate lists
FCcnt = hydro.TcFC(optfolder + "/Tc_subshed"+str(idx))
if landuse == "MRLC Landuse":
STcnt = hydro.TcST_MRLC(optfolder + "/Tc_subshed"+str(idx))
else:
STcnt = hydro.TcST(optfolder + "/Tc_subshed"+str(idx))
if landuse == "MRLC Landuse":
IMPcnt = hydro.TcImpMRLC(optfolder + "/Tc_subshed"+str(idx))
elif landuse == "Ultimate Landuse":
IMPcnt = hydro.TcImpUltimate(optfolder + "/Tc_subshed"+str(idx))
elif landuse == "1997 USGS Landuse" or "1970s USGS Landuse":
IMPcnt = hydro.TcImpUSGS(optfolder + "/Tc_subshed"+str(idx))
else:
IMPcnt = hydro.TcImpAnderson(optfolder + "/Tc_subshed"+str(idx))
FC.append(FCcnt)
ST.append(STcnt)
IA.append(IMPcnt)
lst = lst + 1
# compute percent area of FC, ST, and IA in each subwatershed
FC = [(float(a)/b)*100 for a,b in zip(FC,subarea)]
ST = [(float(a)/b)*100 for a,b in zip(ST,subarea)]
IA = [(float(a)/b)*100 for a,b in zip(IA,subarea)]
# create lists to store gridcode (for duplication of FC, ST, IA, maxlength, and theslope), province, and area
gc_lst = []
prov_lst = []
area_lst = []
sub_prov = arcpy.SearchCursor(subshed_prov, "", "", "GRIDCODE;PROVINCE;Area", "")
for sub in sub_prov:
gc = sub.getValue("GRIDCODE")
gc_lst.append(gc)
prov = sub.getValue("PROVINCE")
prov_lst.append(prov)
area = sub.getValue("Area")
area_lst.append(area/900) # list with subshed area percent in different provinces
del sub_prov
# update length (maxlength) and slope (theslope) fields
longfp = optfolder + "/longfp.dbf"
lfp_max = arcpy.SearchCursor(longfp,"","","MAX","") # already in feets
lfp_list = []
for lfp in lfp_max:
max_lfp = lfp.getValue("MAX")
lfp_list.append(max_lfp)
lfp_lst = 0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"LngFlwPth","FLOAT",15,4)
lfp_update = arcpy.UpdateCursor(subshed_poly)
for lonfp in lfp_update:
sheds_lfp = lfp_list[lfp_lst]
lonfp.LngFlwPth = sheds_lfp
lfp_update.updateRow(lonfp)
lfp_lst = lfp_lst + 1
slp_lst = 0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"Slope","FLOAT","#","#")
slope_update = arcpy.UpdateCursor(subshed_poly)
for slope in slope_update:
sheds_slope = slope_list[slp_lst]
slope.Slope = sheds_slope
slope_update.updateRow(slope)
slp_lst = slp_lst + 1
# adjustment to longest flow path and slope for Tc calculation
maxlength = [float(x)/5280 for x in lfp_list] # maxlength = [float(x)/5280 for x in length]
theslope = [x*5280 for x in slope_list] # theslope = [x*5280 for x in slope]
# update lists to match intersect subshed with poly (for average weighting of Tc)
FC_updated = []
ST_updated = []
IA_updated = []
maxlength_updated = []
theslope_updated = []
tmp_tc=[]
for i, g in enumerate(groupby(gc_lst)):
FC_updated += [FC[i]] * len(list(g[1]))
for i, g in enumerate(groupby(gc_lst)):
ST_updated += [ST[i]] * len(list(g[1]))
for i, g in enumerate(groupby(gc_lst)):
IA_updated += [IA[i]] * len(list(g[1]))
for i, g in enumerate(groupby(gc_lst)):
maxlength_updated += [maxlength[i]] * len(list(g[1]))
for i, g in enumerate(groupby(gc_lst)):
theslope_updated += [theslope[i]] * len(list(g[1]))
# looping over lists to compute "temptc" -- make sure that all lists are of equal length
for a,b,c,d,e,f,g in zip(prov_lst,maxlength_updated,theslope_updated,FC_updated,IA_updated,ST_updated,area_lst):
if a == "A":
temptc = (0.133*(b**(0.475))*(c**(-0.187))*((101-d)**(-0.144))*((101-e)**(0.861))*((f+1)**(0.154))*((10)**(0.194)))
elif a == "W":
temptc = (0.133*(b**(0.475))*(c**(-0.187))*((101-d)**(-0.144))*((101-e)**(0.861))*((f+1)**(0.154))*((10)**(0.366)))
elif a == "E":
temptc = (0.133*(b**(0.475))*(c**(-0.187))*((101-d)**(-0.144))*((101-e)**(0.861))*((f+1)**(0.154))*((10)**(0.366)))
else:
temptc = (0.133*(b**(0.475))*(c**(-0.187))*((101-d)**(-0.144))*((101-e)**(0.861))*((f+1)**(0.154)))
tmp_tc.append(g*temptc)
tc = []
for num, grp in groupby(enumerate(gc_lst), itemgetter(1)):
tmp_list = [tmp_tc[idx] for idx, _ in grp]
tc.append(sum(tmp_list))
tc_list = [a/b for a,b in zip(tc,subarea)]
tc_lst=0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"Tc","FLOAT",15,4)
tc_update = arcpy.UpdateCursor(subshed_poly)
for tc in tc_update:
sheds_tc = tc_list[tc_lst]
tc.Tc = sheds_tc
tc_update.updateRow(tc)
tc_lst = tc_lst + 1
# Velocity Method
else:
# LongPathSub"s schema is locked once added to data frame and "Reset" button fails due to this reason
# If we don"t add LongPathSub to current data frame then we can remove "vel_meth" folder without any error
arcpy.env.addOutputsToMap = False
Tc_method == "Velocity Method Tc Calculation"
#*******************************************************************************************************
# Add fields to "subshed.shp" for subsequent processing during segment merge
#*******************************************************************************************************
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"sheet_n","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"sheet_P","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"sheet_L","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"shal_Paved","TEXT",15,4)
arcpy.AddField_management(subshed_poly,"channel_n","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"ChanDef","TEXT",15,4)
arcpy.AddField_management(subshed_poly,"ChanSA","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"WidthCoef","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"WidthExp","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"DepthCoef","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"DepthExp","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"XAreaCoef","FLOAT",15,4)
arcpy.AddField_management(subshed_poly,"XAreaExp","FLOAT",15,4)
# add "n_sheet"
sheet_n_lst = []
sheet_n_lst.append(Tc_ns)
sheet_n_lst = sheet_n_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for n_sheet in subpoly:
relay = sheet_n_lst[index]
n_sheet.sheet_n = relay
subpoly.updateRow(n_sheet)
index = index + 1
# add "P_sheet"
sheet_P_lst = []
sheet_P_lst.append(Tc_P)
sheet_P_lst = sheet_P_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for P_sheet in subpoly:
relay = sheet_P_lst[index]
P_sheet.sheet_P = relay
subpoly.updateRow(P_sheet)
index = index + 1
# add "L_sheet"
sheet_L_lst = []
sheet_L_lst.append(Tc_L)
sheet_L_lst = sheet_L_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for L_sheet in subpoly:
relay = sheet_L_lst[index]
L_sheet.sheet_L = relay
subpoly.updateRow(L_sheet)
index = index + 1
# add "Paved or Unpaved"
shallow_Paved_lst = []
if Tc_paved == True:
shallow_Paved_lst.append("Paved")
shallow_Paved_lst = shallow_Paved_lst * no_subwatersheds
else:
shallow_Paved_lst.append("Unpaved")
shallow_Paved_lst = shallow_Paved_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for Paved in subpoly:
relay = shallow_Paved_lst[index]
Paved.shal_Paved = relay
subpoly.updateRow(Paved)
index = index + 1
# use "NHD" or modified streams
ChanDef_lst = []
if Tc_NHD == True:
ChanDef_lst.append("NHD")
ChanDef_lst = ChanDef_lst * no_subwatersheds
else:
ChanDef_lst.append("SourceArea")
ChanDef_lst = ChanDef_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for Def_Chan in subpoly:
relay = ChanDef_lst[index]
Def_Chan.ChanDef = relay
subpoly.updateRow(Def_Chan)
index = index + 1
# add "n_channel"
channel_n_lst = []
channel_n_lst.append(Tc_nc)
channel_n_lst = channel_n_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for n_channel in subpoly:
relay = channel_n_lst[index]
n_channel.channel_n = relay
subpoly.updateRow(n_channel)
index = index + 1
# add "SA_sheet"
ChanSA_lst = []
ChanSA_lst.append(Tc_sa)
ChanSA_lst = ChanSA_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for SA_Chan in subpoly:
relay = ChanSA_lst[index]
SA_Chan.ChanSA = relay
subpoly.updateRow(SA_Chan)
index = index + 1
# add "cw_coef"
WidthCoef_lst = []
WidthCoef_lst.append(Tc_cwCoef)
WidthCoef_lst = WidthCoef_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for cw_coef in subpoly:
relay = WidthCoef_lst[index]
cw_coef.WidthCoef = relay
subpoly.updateRow(cw_coef)
index = index + 1
# add "cw_exp"
WidthExp_lst = []
WidthExp_lst.append(Tc_cwExp)
WidthExp_lst = WidthExp_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for cw_exp in subpoly:
relay = WidthExp_lst[index]
cw_exp.WidthExp = relay
subpoly.updateRow(cw_exp)
index = index + 1
# add "cd_coef"
DepthCoef_lst = []
DepthCoef_lst.append(Tc_cdCoef)
DepthCoef_lst = DepthCoef_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for cd_coef in subpoly:
relay = DepthCoef_lst[index]
cd_coef.DepthCoef = relay
subpoly.updateRow(cd_coef)
index = index + 1
# add "cd_exp"
DepthExp_lst = []
DepthExp_lst.append(Tc_cdExp)
DepthExp_lst = DepthExp_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for cd_exp in subpoly:
relay = DepthExp_lst[index]
cd_exp.DepthExp = relay
subpoly.updateRow(cd_exp)
index = index + 1
# add "ca_coef"
XAreaCoef_lst = []
XAreaCoef_lst.append(Tc_caCoef)
XAreaCoef_lst = XAreaCoef_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for ca_coef in subpoly:
relay = XAreaCoef_lst[index]
ca_coef.XAreaCoef = relay
subpoly.updateRow(ca_coef)
index = index + 1
# add "ca_exp"
XAreaExp_lst = []
XAreaExp_lst.append(Tc_caExp)
XAreaExp_lst = XAreaExp_lst * no_subwatersheds
subshed_poly = optfolder + "/subshed.shp"
subpoly = arcpy.UpdateCursor(subshed_poly)
index = 0
for ca_exp in subpoly:
relay = XAreaExp_lst[index]
ca_exp.XAreaExp = relay
subpoly.updateRow(ca_exp)
index = index + 1
# Add fields to "subshed.shp"
longfp = optfolder + "/longfp.dbf"
lfp_max = arcpy.SearchCursor(longfp,"","","MAX","") # already in feets
lfp_list = []
for lfp in lfp_max:
max_lfp = lfp.getValue("MAX")
lfp_list.append(max_lfp)
lfp_lst = 0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"LngFlwPth","FLOAT",15,4)
lfp_update = arcpy.UpdateCursor(subshed_poly)
for lonfp in lfp_update:
sheds_lfp = lfp_list[lfp_lst]
lonfp.LngFlwPth = sheds_lfp
lfp_update.updateRow(lonfp)
lfp_lst = lfp_lst + 1
slope_sheds = optfolder + "/slope_sheds.dbf"
slope_mean = arcpy.SearchCursor(slope_sheds,"","","MEAN","")
slope_list = []
for slp in slope_mean:
mean_slope = slp.getValue("MEAN")
slope_list.append(mean_slope)
slp_lst = 0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"Slope","FLOAT","#","#")
slope_update = arcpy.UpdateCursor(subshed_poly)
for slope in slope_update:
sheds_slope = slope_list[slp_lst]
slope.Slope = sheds_slope
slope_update.updateRow(slope)
slp_lst = slp_lst + 1
# copy "vel_meth" folder from "data" into "optfolder"
src = r"" + Directory + "/data/vel_meth"
dst = optfolder + "/vel_meth"
shutil.copytree(src, dst)
# read sheet global variables [some variables are re-defines just to keep consistency with Avenue code]
sheet_n = float(Tc_ns)
sheet_P = float(Tc_P)
sheet_L = float(Tc_L)
overland_n = sheet_n
overland_P = sheet_P
overland_L = sheet_L
# read channel global variables
chan_n = float(Tc_nc)
thechanSA = float(Tc_sa)
widthcoef = float(Tc_cwCoef)
widthexp = float(Tc_cwExp)
depthcoef = float(Tc_cdCoef)
depthexp = float(Tc_cdExp)
xareacoef = float(Tc_caCoef)
xareaexp = float(Tc_caExp)
# read data
dem = optfolder + "/dem"
dirgrid = optfolder + "/flowdir"
fc = optfolder + "/subshed.shp"
desc = arcpy.Describe(fc)
rows = arcpy.SearchCursor(fc)
# begin indexing and process data for each sub-basin
tc_list = []
for idx,row in enumerate(rows):
aPoly = row.getValue(desc.shapefieldname)
setExtent = aPoly.extent
expPercent = 0.1
bigExtent = hydro.expandExtent(setExtent,expPercent)
arcpy.env.extent = arcpy.Extent(int(bigExtent[0]), int(bigExtent[1]), int(bigExtent[2]), int(bigExtent[3]))
facc = optfolder + "/flowacc"
faccdesc = arcpy.Describe(facc)
cellSize = faccdesc.meanCellHeight
arcpy.env.cellSize = cellSize
# set minimum value
slopetcgrid = optfolder + "/slope_calc/landslope"
slopetcgrid = arcpy.sa.Divide(slopetcgrid, 3.2808) # 07-06-2015: "3.2808" added to factor landslope units into meters
# NOT doing clipping and therefore simply saving same grid with index numbers for time being to avoid
# chaning all of the code below.
# Note: Average slope is now based of whole watershed instead of subwatersheds
slopetcgrid.save(optfolder + "/vel_meth/slopetcgrid" + str(idx)) # 05-09-2017: temporarily added to save grids with index number
dirgrid = optfolder + "/flowdir"
dirgrid = arcpy.sa.Times(dirgrid, 1)
dirgrid.save(optfolder + "/vel_meth/dirgrid" + str(idx))
# create mask raster for each sub-basin
basingrid = optfolder + "/basingrid"
arcpy.MakeFeatureLayer_management(fc, "layer" + str(idx), ' "FID" = ' + str(idx))
mask = arcpy.Clip_management(basingrid, "#", optfolder + "/vel_meth/mask" + str(idx), "layer" + str(idx), "0", "ClippingGeometry")
# clip and save sub-sheds 10% bigger than the original extent
avgslope_val = arcpy.GetRasterProperties_management(optfolder + "/vel_meth/slopetcgrid" + str(idx),"MEAN")
avgslope = float(avgslope_val.getOutput(0))
# mask above clipped layers -- make sure to only mask those at this point which will not
# be used in processing below (only mask which will be converted to ascii straight away)
facc = arcpy.sa.Plus(facc,1)
facc.save(optfolder + "/vel_meth/areagrid" + str(idx))
elev = arcpy.sa.Times(dem,mask)
elev.save(optfolder + "/vel_meth/elevgrid" + str(idx))
if Tc_NHD == True:
# 05-09-2017: temporarily disabling clipping of NHD streams
arcpy.PolylineToRaster_conversion(Directory + "/data/maryland/nhd_streamsm.shp", "FID", optfolder + "/vel_meth/streamgrid" + str(idx), "MAXIMUM_LENGTH", "NONE", "30")
strmgr_flt = arcpy.sa.Float(optfolder + "/vel_meth/streamgrid" + str(idx))
strmgr_flt.save(optfolder + "/vel_meth/strmgr_flt" + str(idx))
if Tc_infStreams == True:
srcpixel = (thechanSA * pow(5280,2))/pow((3.28084 * cellSize),2)
streamgrid_con = arcpy.sa.Con(facc >= srcpixel, 1, 0)
InfStr_null = arcpy.sa.SetNull(streamgrid_con, streamgrid_con, "Value = 0")
InfStr_min = arcpy.sa.Minus(InfStr_null,1)
streamgrid_masked = arcpy.sa.Times(InfStr_min,mask)
streamgrid_masked.save(optfolder + "/vel_meth/streamgrid" + str(idx))
strmgr_flt = arcpy.sa.Float(optfolder + "/vel_meth/streamgrid" + str(idx))
strmgr_flt.save(optfolder + "/vel_meth/strmgr_flt" + str(idx))
# get upstream and log2 base direction raster for sub-extents
upgrid = arcpy.sa.FlowLength(optfolder + "/vel_meth/dirgrid" + str(idx), "UPSTREAM", "")
upgrid = arcpy.sa.Times(upgrid, 3.28084)
upgrid = arcpy.sa.Times(upgrid,mask)
upgrid.save(optfolder + "/vel_meth/upgrid" + str(idx))
dlgrid_temp1 = arcpy.sa.Log2(optfolder + "/vel_meth/dirgrid" + str(idx))
dlgrid_temp2 = dlgrid_temp1 % 2
dlgrid_temp3 = arcpy.sa.Con(dlgrid_temp2 > 0, pow(2,0.5), 1)
dlgrid_temp4 = arcpy.sa.Times(dlgrid_temp3,cellSize)
dlgrid = arcpy.sa.Times(dlgrid_temp4, 3.28084)
dlgrid = arcpy.sa.Times(dlgrid,mask)
dlgrid.save(optfolder + "/vel_meth/dlgrid" + str(idx))
upgridp = upgrid + dlgrid
upgridp.save(optfolder + "/vel_meth/upgridp" + str(idx))
# calculate indicator grid ("indic = 3" means everything is swale)
indic = arcpy.sa.Times(basingrid,3) # create swale raster
indic = arcpy.sa.Con(upgridp <= sheet_L, 1, indic) # substitute sheet flow
indic = arcpy.sa.Con((upgridp > sheet_L) & (upgrid < sheet_L), 2, indic) # substitute pixels that are part sheet, part swale
channel_con = arcpy.sa.IsNull(optfolder + "/vel_meth/streamgrid" + str(idx))
indic = arcpy.sa.Con(channel_con == 0, 4, indic) # substitute channel pixels
indic = arcpy.sa.Times(indic,mask)
indic.save(optfolder + "/vel_meth/indic" + str(idx))
# 5/8/2017: condition to check if slope is <= 0 then use value 0.001 which is roughly approximates a 0.1 foot drop over a 30 meter pixel
slopegrid = arcpy.sa.Con(optfolder + "/vel_meth/slopetcgrid" + str(idx) <= 0, 0.001, optfolder + "/vel_meth/slopetcgrid" + str(idx))
slopegrid = arcpy.sa.IsNull(slopegrid)
slopegrid = arcpy.sa.Con(slopegrid == 1, 1/(1.4142 * cellSize * 3.28084), optfolder + "/vel_meth/slopetcgrid" + str(idx))
slopegrid = arcpy.sa.Con((indic <= 2) & (slopegrid == 0), avgslope, slopegrid) # ^ this has already been taken care off
slopegrid = arcpy.sa.Con(slopegrid <= 0, 0.001, slopegrid)
slopegrid = arcpy.sa.Times(slopegrid,mask)
slopegrid.save(optfolder + "/vel_meth/slopegrid" + str(idx))
ol_temp1 = arcpy.sa.Minus(sheet_L, upgrid)
ol_temp2 = arcpy.sa.Con(indic == 2, ol_temp1, 0)
ol = arcpy.sa.Divide(arcpy.sa.Power(ol_temp2, 0.8), arcpy.sa.Power(dlgrid, 0.8))
sl_temp1 = arcpy.sa.Minus(upgridp, sheet_L)
sl_temp2 = arcpy.sa.Con(indic == 2, sl_temp1, 0)
sl = arcpy.sa.Divide(sl_temp2, dlgrid)
# calcualte Sheet part of travel time
sheet_Tt_tm1 = arcpy.sa.Times(sheet_n, dlgrid)
sheet_Tt_tm2 = pow(sheet_P, 0.5) # Since no raster involved so should not use "arcpy.sa.Power"
sheet_Tt_tm3 = arcpy.sa.Power(slopegrid, 0.4)
sheet_Tt_tm4 = arcpy.sa.Power(sheet_Tt_tm1, 0.8)
sheet_Tt_tm5 = arcpy.sa.Times(0.007, sheet_Tt_tm4)
sheet_Tt_tm6 = arcpy.sa.Times(sheet_Tt_tm2, sheet_Tt_tm3)
sheet_Tt = arcpy.sa.Divide(sheet_Tt_tm5, sheet_Tt_tm6) # 07-16-2015: Eq. 3-3 in TR-55 Document
tt_sheet = arcpy.sa.Con(indic == 1, sheet_Tt, 0)
tt_sheet.save(optfolder + "/vel_meth/tt_sheet" + str(idx))
# calculate Swale part of travel time
if Tc_paved:
swale_coef = 73181.52 # revised from 73440 based on 3600 times the value in Appendix F of TR-55 Manual
if Tc_unpaved:
swale_coef = 58084.20 # revised from 57600 based on 3600 times the value in Appendix F of TR-55 Manual
slopegrid_sqrt = arcpy.sa.SquareRoot(optfolder + "/vel_meth/slopegrid" + str(idx))
swale_Tt = arcpy.sa.Divide(dlgrid, (slopegrid_sqrt * swale_coef))
tt_swale = arcpy.sa.Con(indic == 3, swale_Tt, 0) # 5/7/2017: separating sheet, swale, mixed, and channel rasters based on indicator grid
tt_swale.save(optfolder + "/vel_meth/tt_swale" + str(idx))
mixed_tt_temp1 = arcpy.sa.Times(ol,sheet_Tt)
mixed_tt_temp2 = arcpy.sa.Times(sl,swale_Tt)
mixed_tt = arcpy.sa.Plus(mixed_tt_temp1, mixed_tt_temp2)
tt_mixed = arcpy.sa.Con(indic == 2, mixed_tt, 0) # 5/7/2017: separating sheet, swale, mixed, and channel rasters based on indicator grid
tt_mixed.save(optfolder + "/vel_meth/tt_mixed" + str(idx))
const = (pow(cellSize,2) / 27878400) * (pow(3.28084,2)) # 07-09-2015: Added brackets to divide only with 27878400
areagrid_const = arcpy.sa.Times(optfolder + "/vel_meth/areagrid" + str(idx),const)
areagrid_xarea = arcpy.sa.Power(areagrid_const,xareaexp)
chanarea = arcpy.sa.Times(xareacoef, areagrid_xarea)
chanarea.save(optfolder + "/vel_meth/chanarea" + str(idx))
areagrid_width = arcpy.sa.Power(areagrid_const,widthexp)
chanwidth = arcpy.sa.Times(widthcoef, areagrid_width)
chanwidth.save(optfolder + "/vel_meth/chanwidth" + str(idx))
areagrid_depth = arcpy.sa.Power(areagrid_const,depthexp)
chandepth = arcpy.sa.Times(depthcoef, areagrid_depth)
chandepth.save(optfolder + "/vel_meth/chandepth" + str(idx))
chanperim = arcpy.sa.Plus((arcpy.sa.Times(2, chandepth)), chanwidth)
Rh = arcpy.sa.Divide(chanarea, chanperim)
chan_vel_temp1 = arcpy.sa.Power(Rh,0.666667)
chan_vel_1 = arcpy.sa.Divide(1.49, chan_n)
chan_vel_2 = arcpy.sa.Times(chan_vel_temp1, slopegrid_sqrt)
chan_vel_3 = arcpy.sa.Times(chan_vel_1, chan_vel_2)
chan_vel = arcpy.sa.Times(chan_vel_3, 3600) # 07-07-2015: Eq. 3-4 in TR-55 Document
chan_tt = arcpy.sa.Divide(dlgrid, chan_vel) # "chan_Tt" is same as "chan_tt" [Avenue is case insensitive]
chan_tt.save(optfolder + "/vel_meth/chan_tt" + str(idx))
# new chan_tt grid created as per new cumulatettnew.exe requirements
chan_tt_new = arcpy.sa.Con(indic == 4, 0.001, 0)
chan_tt_new.save(optfolder + "/vel_meth/chan_tt_new" + str(idx))
# add all tt rasters
tt_new1 = arcpy.sa.Plus(optfolder + "/vel_meth/tt_swale" + str(idx),optfolder + "/vel_meth/tt_mixed" + str(idx))
tt_new2 = arcpy.sa.Plus(optfolder + "/vel_meth/tt_swale" + str(idx),optfolder + "/vel_meth/chan_tt_new" + str(idx))
tt_new3 = arcpy.sa.Plus(tt_new1,tt_new2)
tt_new3.save(optfolder + "/vel_meth/tt_new" + str(idx))
strm_flt = arcpy.sa.IsNull(optfolder + "/vel_meth/strmgr_flt" + str(idx))
tt = arcpy.sa.Con(strm_flt == 0, chan_tt, tt_mixed)
tt.save(optfolder + "/vel_meth/tt" + str(idx))
dirgrid_moglen = arcpy.sa.Reclassify(optfolder + "/vel_meth/dirgrid" + str(idx),"Value", "1 1;2 8;4 7;8 6;16 5;32 4;64 3;128 2","DATA")
dirgrid_moglen.save(optfolder + "/vel_meth/dirgrid_mog" + str(idx)) # 5/7/2017: using big extent
# create slopetc grid masked files -- need original "slope_calc" file as input for cumulatettnew.exe
slopetc_grid = optfolder + "/vel_meth/slopetcgrid" + str(idx)
slopetc_grid = arcpy.sa.Times(slopetc_grid,mask)
slopetc_grid.save(optfolder + "/vel_meth/slptcgrid" + str(idx))
areagrid_mask = arcpy.sa.Times(optfolder + "/vel_meth/areagrid" + str(idx),mask)
areagrid_mask.save(optfolder + "/vel_meth/areagr_new" + str(idx))
chanarea_new = arcpy.sa.Times(optfolder + "/vel_meth/chanarea" + str(idx),mask)
chanarea_new.save(optfolder + "/vel_meth/chanarea_n" + str(idx))
chanwidth_new = arcpy.sa.Times(optfolder + "/vel_meth/chanwidth" + str(idx),mask)
chanwidth_new.save(optfolder + "/vel_meth/chanwidth_n" + str(idx))
chandepth_new = arcpy.sa.Times(optfolder + "/vel_meth/chandepth" + str(idx),mask)
chandepth_new.save(optfolder + "/vel_meth/chandepth_n" + str(idx))
arcpy.env.extent = optfolder + "/vel_meth/elevgrid" + str(idx)
# create ascii files
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/dirgrid_mog" + str(idx), optfolder + "/vel_meth/temp1.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/areagr_new" + str(idx), optfolder + "/vel_meth/temp2.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/tt_new" + str(idx), optfolder + "/vel_meth/temp3.asc") # 5/7/2017: new tt
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/indic" + str(idx), optfolder + "/vel_meth/temp6.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/chanwidth_n" + str(idx), optfolder + "/vel_meth/temp7.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/chandepth_n" + str(idx), optfolder + "/vel_meth/temp8.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/chanarea_n" + str(idx), optfolder + "/vel_meth/temp9.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/dlgrid" + str(idx), optfolder + "/vel_meth/temp10.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/slptcgrid" + str(idx), optfolder + "/vel_meth/temp11.asc")
arcpy.RasterToASCII_conversion(optfolder + "/vel_meth/elevgrid" + str(idx), optfolder + "/vel_meth/temp12.asc")
# write string variables for "cumulate.txt"
filestring = ""
filestring = filestring + "temp1.asc" + "\n"
filestring = filestring + "temp2.asc" + "\n"
filestring = filestring + "temp3.asc" + "\n"
filestring = filestring + "temp6.asc" + "\n"
filestring = filestring + "temp7.asc" + "\n"
filestring = filestring + "temp8.asc" + "\n"
filestring = filestring + "temp9.asc" + "\n"
filestring = filestring + "temp10.asc" + "\n"
filestring = filestring + "temp11.asc" + "\n"
filestring = filestring + "temp12.asc" + "\n"
filestring = filestring + "temp5.asc" + "\n"
filestring = filestring + "temp4.asc" + "\n"
filestring = filestring + "velmethtable" + str(idx) + ".txt" + "\n"
filestring = filestring + str(overland_L) + "\n"
filestring = filestring + str(overland_n) + "\n"
filestring = filestring + str(overland_P) + "\n"
filestring = filestring + str(swale_coef) + "\n"
filestring = filestring + str(chan_n)
# write above order to a text file -- "cumulate.txt"
defFN = optfolder + "/vel_meth/cumulate.txt"
cumulate = open(defFN, "w")
cumulate.write(filestring)
cumulate.close()
# run "cumulatettnew.exe" to generate "velmethtables.txt"
os.chdir(optfolder + "/vel_meth")
os.system(optfolder + "/vel_meth/cumulatettnew.exe")
time.sleep(4)
# convert "longfile" [temp5.asc] to GRID, generate its attribute table and join with velmeth tables
arcpy.ASCIIToRaster_conversion(optfolder + "/vel_meth/temp5.asc", optfolder + "/vel_meth/LongPathSub" + str(idx), "INTEGER")
arcpy.BuildRasterAttributeTable_management(optfolder + "/vel_meth/LongPathSub" + str(idx))
# Delete ASCII files to run executable in next iteration of for loop without error
arcpy.Delete_management(optfolder + "/vel_meth/temp1.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp2.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp3.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp4.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp5.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp6.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp7.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp8.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp9.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp10.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp11.asc")
arcpy.Delete_management(optfolder + "/vel_meth/temp12.asc")
# create "tc_list" ["tc_list" is defined at the start of this loop]
velmethtables = optfolder + "/vel_meth/velmethtable"
infile = open(velmethtables + str(idx) + ".txt","r").readlines()
theline = infile[-1].split()
Tot_Time = theline[13]
tc_list.append(Tot_Time)
del idx
del row
# add "tc_list" to "subshed.shp"
tc_lst=0
subshed_poly = optfolder + "/subshed.shp"
arcpy.AddField_management(subshed_poly,"Tc","FLOAT",15,4)
tc_update = arcpy.UpdateCursor(subshed_poly)
for tc in tc_update:
sheds_tc_list = tc_list[tc_lst]
sheds_tc = float(sheds_tc_list) # 5/23/2017: type is adjusted to match above Field type
tc.Tc = sheds_tc
tc_update.updateRow(tc)
tc_lst = tc_lst + 1
# join "longfile" [temp5.asc] with velmeth tables
fc = optfolder + "/subshed.shp"
desc = arcpy.Describe(fc)
subs = arcpy.SearchCursor(fc)
arcpy.env.qualifiedFieldNames = False # to have attribute tables with original names
fields = ["Pixel", "Type", "Mixed", "Drain_Area", "Elev", "Slope", "Width", "Depth", "Xarea",
"I_Length", "Tot_Length", "Vel.", "I_Time", "Tot_Time"]
global arcid_list
global arcid
subshed = optfolder + "/subshed.shp"
arcid_list = []
shedtab = arcpy.SearchCursor(subshed,"","","ARCID","")
for s in shedtab:
arcid = s.getValue("ARCID")
arcid_list.append(str(int(arcid)))
arcpy.env.addOutputsToMap = False
arcpy.env.extent = "MAXOF"
for idx,i in enumerate(subs):
arcpy.TableToDBASE_conversion(optfolder + "/vel_meth/velmethtable" + str(idx) + ".txt", optfolder + "/vel_meth")
arcpy.JoinField_management(optfolder + "/vel_meth/LongPathSub" + str(idx), "value", optfolder + "/vel_meth/velmethtable" + str(idx) + ".dbf", "Pixel", fields)
del fc, desc, subs, fields, idx, i
# re-reading because above variables lead to shapefile lock and prohibit file deletion
fc = optfolder + "/subshed.shp"
subs = arcpy.SearchCursor(fc)
arcpy.env.qualifiedFieldNames = False # to have attribute tables with original names
fields = ["Pixel", "Type", "Mixed", "Drain_Area", "Elev", "Slope", "Width", "Depth", "Xarea",
"I_Length", "Tot_Length", "Vel.", "I_Time", "Tot_Time"]
for s,sub in enumerate(subs):
arcpy.env.addOutputsToMap = False
lp_rc = arcpy.sa.Reclassify(optfolder + "/vel_meth/LongPathSub" + str(s), "TYPE", RemapValue([["",0],["channel",1],["swale",2],["overland",3],["mixed",4]]),"DATA")
lp_rc.save(optfolder + "/vel_meth/lp" + str(s) + "_reclass")
lp_sn = arcpy.sa.SetNull(optfolder + "/vel_meth/lp" + str(s) + "_reclass", optfolder + "/vel_meth/lp" + str(s) + "_reclass", "Value = 0")
lp_sn.save(optfolder + "/vel_meth/lp" + str(s) + "_setnull")
arcpy.env.addOutputsToMap = True
arcpy.RasterToPolyline_conversion(optfolder + "/vel_meth/lp" + str(s) + "_setnull", optfolder + "/aux_folder/Longest_Path_Sub_" + str(s+1) + ".shp","ZERO","0","NO_SIMPLIFY","Value")
del fc, subs, s, lp_rc, lp_sn
#*******************************************************************************************************
# turn layers ON/OFF in current data frame
#*******************************************************************************************************
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
layers = arcpy.mapping.ListLayers(mxd, "", df)
for lyr in layers:
if lyr.name == "subshed_prov":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lp": # from above step we likely get a single layer in memory over and over so have to remove it
arcpy.mapping.RemoveLayer(df, lyr)
for i in xrange(0,no_subwatersheds,1):
if lyr.name == "Tc_subshed"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "Tc_temp"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "shedMask"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "layer"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "streamgrid"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "dirgrid"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "mask"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "slopetcgrid"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "nhd_str"+str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lpsub" + str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lp" + str(i) + "_reclass":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "lp" + str(i) + "_setnull":
arcpy.mapping.RemoveLayer(df, lyr)
if Tc_method == "Velocity Method Tc Calculation":
if lyr.name == "LongPathSub" + str(i):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "Longest_Path_Sub_" + str(i+1) + "":
lyr.visible = True
arcpy.ApplySymbologyFromLayer_management(lyr,r"" + Directory + "/data/mdfiles/legends/longpathsub.lyr")
if lyr.name == "elevmerge":
lyr.visible = False
if lyr.name == "elevzones":
lyr.visible = False
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
del lyr
del i
#*******************************************************************************************************
# clean "vel_meth" directory after indexed looping [delete everything except "velmethtables" and *.exe]
#*******************************************************************************************************
for i in xrange(0,no_subwatersheds,1):
if os.path.exists(optfolder + "/vel_meth/areagrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/areagrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/areagr_new" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/areagr_new" + str(i))
if os.path.exists(optfolder + "/vel_meth/chanarea_n" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/chanarea_n" + str(i))
if os.path.exists(optfolder + "/vel_meth/chanarea" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/chanarea" + str(i))
if os.path.exists(optfolder + "/vel_meth/chandepth_n" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/chandepth_n" + str(i))
if os.path.exists(optfolder + "/vel_meth/chandepth" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/chandepth" + str(i))
if os.path.exists(optfolder + "/vel_meth/chanwidth_n" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/chanwidth_n" + str(i))
if os.path.exists(optfolder + "/vel_meth/chanwidth" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/chanwidth" + str(i))
if os.path.exists(optfolder + "/vel_meth/dirgrid_mog" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/dirgrid_mog" + str(i))
if os.path.exists(optfolder + "/vel_meth/dirgrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/dirgrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/dlgrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/dlgrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/elevgrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/elevgrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/indic" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/indic" + str(i))
if os.path.exists(optfolder + "/vel_meth/mask" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/mask" + str(i))
if os.path.exists(optfolder + "/vel_meth/slopegrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/slopegrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/slopetcgrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/slopetcgrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/slptcgrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/slptcgrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/streamgrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/streamgrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/tt" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/tt" + str(i))
if os.path.exists(optfolder + "/vel_meth/upgrid" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/upgrid" + str(i))
if os.path.exists(optfolder + "/vel_meth/upgridp" + str(i)):
arcpy.Delete_management(optfolder + "/vel_meth/upgridp" + str(i))
if os.path.exists(optfolder + "/vel_meth/lp" + str(i) + "_setnull"):
arcpy.Delete_management(optfolder + "/vel_meth/lp" + str(i) + "_setnull")
if os.path.exists(optfolder + "/vel_meth/lp" + str(i) + "_reclass"):
arcpy.Delete_management(optfolder + "/vel_meth/lp" + str(i) + "_reclass")
pythonaddins.MessageBox("Attributes Calculation Complete!", "Calculate Attributes", 0)
#*******************************************************************************************************
# turn subwatershed delineation OFF and Stream cross-section or Precipitation Depths tool ON
#*******************************************************************************************************
button10.enabled = False
FromNode = []
ToNode = []
subriver = optfolder + "/subrivers.shp"
sr = arcpy.SearchCursor(subriver, "", "", "GRID_CODE;FROM_NODE;TO_NODE", "")
for node in sr:
fn = node.getValue("FROM_NODE")
tn = node.getValue("TO_NODE")
FromNode.append(fn)
ToNode.append(tn)
# "reach_check_lst" is to check if there any reaches.
# If there are no reaches then watershed is treated as with single basin.
reach_check_lst = [x for x in FromNode if x in ToNode]
if not reach_check_lst:
tool6.enabled = False # changed from tool.5 to tool.8
button13.enabled = True
else:
tool6.enabled = True # changed from tool.5 to tool.8
button13.enabled = False
button11.enabled = True
button12.enabled = True
# If "Velocity method" is used then we need combine longest flow path,
# otherwise grey out button
if Tc_method == "Velocity Method Tc Calculation":
button11.enabled = True
else:
button11.enabled = False
arcpy.env.extent = "MAXOF"