GISHydroNXT System Documentation
#*********************************************************************************
# Author: UMD
# Date: 24-07-2018
# Modified: n/a
# Classes: TransectLine() ; Xeditor()
# Functions: n/a
# Modules: arcpy
# Comments: n/a
#*********************************************************************************
class TransectLine(object):
"""Implementation for GISHydroNXT_addin.tool6 (Tool)"""
def __init__(self):
self.enabled = False
self.cursor = 1
self.shape = "Line"
def onLine(self, line_geometry):
arcpy.env.scratchWorkspace = scratchfolder
arcpy.env.workspace = optfolder
arcpy.env.extent = "MAXOF"
#*******************************************************************************************************
# Add DEM elevation and distance fields to transect line
#*******************************************************************************************************
arcpy.CopyFeatures_management(line_geometry, optfolder + "/trans_line.shp") # saving transect line
with arcpy.da.SearchCursor(optfolder + "/trans_line.shp", "SHAPE@") as rows:
for row in rows:
if row[0].length > 500:
pythonaddins.MessageBox("Transect length is "+str(round((row[0].length),1))+" meters. Maximum allowable length is 500m. Please draw it again.","Transect Length Warning",0)
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
layers = arcpy.mapping.ListLayers(mxd, "", df)
for lyr in layers:
if lyr.name == "trans_line":
arcpy.mapping.RemoveLayer(df, lyr)
arcpy.Delete_management(optfolder + "/trans_line.shp","")
break
else:
line_trun = [optfolder + "/trans_line.shp", optfolder + "/subshed.shp"]
arcpy.Intersect_analysis(line_trun, optfolder + "/line_trun.shp", "ALL","#","INPUT")
arcpy.Clip_management(optfolder + "/flowacc", "#", optfolder + "/flowacc_x", optfolder + "/line_trun.shp", "ClippingGeometry")
arcpy.Clip_management(optfolder + "/dem", "#", optfolder + "/dem_x", optfolder + "/line_trun.shp", "ClippingGeometry")
arcpy.CreateRandomPoints_management(optfolder, "transect", optfolder + "/line_trun.shp", "#", "#", 6, "#", "#") # 02-15-2015: changed sampling interval from 3 to 6
arcpy.AddField_management(optfolder + "/transect.shp","Distance","FLOAT",15,4)
# add distance along transect
rows = arcpy.UpdateCursor(optfolder + "/transect.shp")
distance = 0
count=0
for row in rows:
count=count+1
distance = distance+9.8424
row.Distance = distance
rows.updateRow(row)
# delete field "CID"
arcpy.DeleteField_management(optfolder + "/transect.shp", "CID")
# extract multivalues to points using DEM
arcpy.sa.ExtractMultiValuesToPoints(optfolder + "/transect.shp", optfolder + "/dem_x", "NONE")
arcpy.sa.ExtractMultiValuesToPoints(optfolder + "/transect.shp", optfolder + "/flowacc_x", "NONE")
# correct start to end distance and elevation values in feet
rows = arcpy.UpdateCursor(optfolder + "/transect.shp")
for row in rows:
row.setValue("Distance", row.getValue("Distance") - 9.8424)
row.setValue("dem_x", row.getValue("dem_x")) # elevation values are already in feet -- no need for conversion
rows.updateRow(row)
dem_list = []
distance_list = []
flowacc_list = []
upstreamDA = optfolder + "/transect.shp"
usda = arcpy.SearchCursor(upstreamDA, "", "", "dem_x;Distance;flowacc_x", "")
for da in usda:
dem = da.getValue("dem_x")
dem_list.append(dem)
distance = da.getValue("Distance")
distance_list.append(distance)
flow = da.getValue("flowacc_x")
flowacc_list.append(flow)
global trans_distance
trans_distance = [round(elem,1) for elem in distance_list]
global trans_dem
trans_dem = [round(elem,1) for elem in dem_list]
#*******************************************************************************************************
# Transection line width, min, and max elevation
# TWL = Transection Width Length
# TWE_min = minimum Transection Width Elevation
# TWE_max = maximum Transection Width Elevation
#*******************************************************************************************************
global TWL
TWL_max = float(max(distance_list))
TWL = "{0:.2f}".format(TWL_max)
global TWE_min
TWE1 = float(min(dem_list))
TWE_min = "{0:.2f}".format(TWE1)
global TWE_max
TWE2 = float(max(dem_list))
TWE_max = "{0:.2f}".format(TWE2)
#*******************************************************************************************************
# Calculate upstream drainage area
#*******************************************************************************************************
upda_max = max(flowacc_list)
flowacc_cell = optfolder + "/flowacc"
rast = arcpy.Raster(flowacc_cell)
cellsize = rast.meanCellWidth
cellsq = cellsize* cellsize
global areami2_usda
areami2_usda = float((upda_max* cellsq) / 2588881) # conversion into sq miles
areami2_usda = "{0:.2f}".format(areami2_usda)
#*******************************************************************************************************
# get reach slope using subsheds, transect, and subriver shapefiles
#*******************************************************************************************************
arcpy.env.qualifiedFieldNames = False # to have attribute tables with original names
rSlope_intersect = [optfolder + "/line_trun.shp", optfolder + "/subshed.shp"]
arcpy.Intersect_analysis(rSlope_intersect, optfolder + "/reachslope.shp", "ALL","#","INPUT")
arcpy.DeleteField_management(optfolder + "/reachslope.shp", "FID_trans_;Id;AreaMi2;TcMethod;CurveNum;LngFlwPth;Tc;Slope")
arcpy.MakeFeatureLayer_management(optfolder + "/reachslope.shp", "reachlayer")
FID = "ARCID"
Code = "ARCID"
arcpy.AddJoin_management("reachlayer", FID, optfolder + "/subrivers.shp", Code)
arcpy.FeatureClassToShapefile_conversion("reachlayer", optfolder)
reach = arcpy.SearchCursor(optfolder + "/reachlayer.shp", "", "", "Slope", "")
for row in reach:
reach_slope = row.getValue("Slope")
del row
del reach
global reachslope
reachslope = "{0:.5f}".format(reach_slope)
#*******************************************************************************************************
# get bankfull channel width and depth
#*******************************************************************************************************
sumarea = 0
theVTab = optfolder + "/theVTab.dbf"
theVTab = arcpy.SearchCursor("theVTab","","","Count","")
for each in theVTab:
count = each.getValue("Count")
sumarea = sumarea+count
sumArea = sumarea
del each
AParea = 0
PDarea = 0
CParea = 0
theVTab = optfolder + "/theVTab.dbf"
theVTab = arcpy.SearchCursor("theVTab","","","Province;Count","")
for row in theVTab:
theArea = float(row.getValue("Count"))
areapercent = float((theArea/sumArea)* 100)
if row.getValue("Province") == "A":
AParea = AParea + float(areapercent)
elif row.getValue("Province") == "B":
AParea = AParea + float(areapercent)
elif row.getValue("Province") == "P":
PDarea = PDarea + float(areapercent)
elif row.getValue("Province") == "W":
CParea = CParea + float(areapercent)
elif row.getValue("Province") == "E":
CParea = CParea + float(areapercent)
del row
regionarea = []
AParea = AParea/100
PDarea = PDarea/100
CParea = CParea/100
regionarea.append(AParea)
regionarea.append(PDarea)
regionarea.append(CParea)
a = [13.87, 14.78, 10.3]
b = [0.44, 0.39, 0.38]
c = [0.95, 1.18, 1.01]
d = [0.31, 0.34, 0.32]
Coef_W = sum(map(mul, regionarea, a))
Exp_W = sum(map(mul, regionarea, b))
global Wbf
WidthBF = Coef_W* (float(areami2_usda)** Exp_W)
Wbf = "{0:.2f}".format(WidthBF)
Coef_D = sum(map(mul, regionarea, c))
Exp_D = sum(map(mul, regionarea, d))
global Dbf
DepthBF = Coef_D* (float(areami2_usda)** Exp_D)
Dbf = "{0:.2f}".format(DepthBF)
#*******************************************************************************************************
# 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 == "transect":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "reachlayer":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "reachslope":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "elevzones":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "dem_x":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "flowacc_x":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "elevmerge":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "line_trun":
arcpy.mapping.RemoveLayer(df, lyr)
if os.path.exists(optfolder + "/reachlayer.shp"):
arcpy.Delete_management(optfolder + "/reachlayer.shp","")
#*******************************************************************************************************
# call Cross-Section Editor dialog
#*******************************************************************************************************
dlg4 = Xeditor()
class Xeditor(wx.Frame):
def __init__(self):
wx.Frame.__init__(self, None, -1, "Cross Section Editor", size=(590, 390))
self.Bind(wx.EVT_CLOSE, self.OnClose)
panel = wx.Panel(self, -1)
wx.StaticBox(panel, -1, "Transect Line Geometry", (20, 15), size=(260,120))
wx.StaticText(panel, -1, "Transect Line Width:", (40, 35))
wx.StaticText(panel, -1, str(TWL), (200, 35))
wx.StaticText(panel, -1, "ft", (255, 35))
wx.StaticText(panel, -1, "Maximum Elevation:", (40, 60))
wx.StaticText(panel, -1, str(TWE_max), (200, 60))
wx.StaticText(panel, -1, "ft", (255, 60))
wx.StaticText(panel, -1, "Minimum Elevation:", (40, 85))
wx.StaticText(panel, -1, str(TWE_min), (200, 85))
wx.StaticText(panel, -1, "ft", (255, 85))
wx.StaticText(panel, -1, "Upstream Discharge Area:", (40, 110))
wx.StaticText(panel, -1, str(areami2_usda), (200, 110))
wx.StaticText(panel, -1, "mi^2", (245, 110))
wx.StaticBox(panel, -1, "Reach Characteristics", (20, 140), size=(260,90))
wx.StaticText(panel, -1, "Reach Slope", (40, 170))
self.reachslope = wx.TextCtrl(panel, -1, value=str(reachslope), pos=(150, 165), size=(75,25))
wx.StaticText(panel, -1, "ft/ft", (245, 170))
wx.StaticText(panel, -1, "Bankfull Elevation", (40, 200))
self.Ebf = wx.TextCtrl(panel, -1, value=str(TWE_min), pos=(150, 195), size=(75,25)) # Ebf = Bank Full Elevation
wx.StaticText(panel, -1, "ft", (245, 200))
wx.StaticBox(panel, -1, "Channel Geometry", (20, 235), size=(260,90))
wx.StaticText(panel, -1, "Bankful Channel Width", (40, 260))
self.Wbf = wx.TextCtrl(panel, -1, value=str(Wbf), pos=(180, 255), size=(70,25)) # Wbf = Bank Full Channel Width
wx.StaticText(panel, -1, "ft", (260, 260))
wx.StaticText(panel, -1, "Bankful Channel Depth", (40, 290))
self.Dbf = wx.TextCtrl(panel, -1, value=str(Dbf), pos=(180, 285), size=(70,25)) # Dbf = Bank Full Channel Depth
wx.StaticText(panel, -1, "ft", (260, 290))
wx.StaticBox(panel, -1, "Roughness Characteristics", (290, 15), size=(260,140))
wx.StaticText(panel, -1, "Main Channel n Value", (310, 35))
self.nMainValue = wx.TextCtrl(panel, -1, value="0.050", pos=(460, 30), size=(70,25)) # Main Channel n Value
wx.StaticText(panel, -1, "Left* Overbank n Value", (310, 65))
self.nLeftValue = wx.TextCtrl(panel, -1, value="0.100", pos=(460, 60), size=(70,25)) # Left Overbank n Value
wx.StaticText(panel, -1, "Right* Overbank n Value", (310, 95))
self.nRightValue = wx.TextCtrl(panel, -1, value="0.100", pos=(460, 90), size=(70,25)) # Right Overbank n Value
wx.StaticText(panel, -1, "* Facing Downstream", (400, 120))
self.rb1 = wx.RadioButton(panel, -1, "Calculate from GIS Data", (390,175), style=wx.RB_GROUP)
self.rb2 = wx.RadioButton(panel, -1, "Load rating table from file", (390,200))
self.btnRecalculate = wx.Button(panel, label="Recalculate", pos=(290,180))
self.Bind(wx.EVT_BUTTON, self.OnRecalculate, id = self.btnRecalculate.GetId())
self.btnRecalculate.Disable() # 1/3/2017: Activate only if user try to type in any of options in dialog box
self.btnExport = wx.Button(panel, label="Export Cross Section", pos=(290,240))
self.Bind(wx.EVT_BUTTON, self.OnExport, id = self.btnExport.GetId())
self.btnPlot = wx.Button(panel, label="Plot Cross Section", pos=(440,240))
self.Bind(wx.EVT_BUTTON, self.OnPlot, id = self.btnPlot.GetId())
self.btnApply = wx.Button(panel, label="Apply", pos=(290,290))
self.Bind(wx.EVT_BUTTON, self.OnApply, id = self.btnApply.GetId())
self.btnClose = wx.Button(panel, label="Close", pos=(440,290))
self.Bind(wx.EVT_BUTTON, self.OnClose, id = self.btnClose.GetId())
self.Show(True)
self.Centre(True)
style = self.GetWindowStyle()
self.SetWindowStyle( style | wx.STAY_ON_TOP )
self.SetWindowStyle( style )
def OnClose(self, event):
#******************************************************************************************************
# delete layers to avoid duplication conflict
#******************************************************************************************************
if os.path.exists(optfolder + "/trans_line.shp"):
arcpy.Delete_management(optfolder + "/trans_line.shp","")
if os.path.exists(optfolder + "/flowacc_x"):
arcpy.Delete_management(optfolder + "/flowacc_x","")
if os.path.exists(optfolder + "/dem_x"):
arcpy.Delete_management(optfolder + "/dem_x","")
if os.path.exists(optfolder + "/transect.shp"):
arcpy.Delete_management(optfolder + "/transect.shp","")
if os.path.exists(optfolder + "/reachslope.shp"):
arcpy.Delete_management(optfolder + "/reachslope.shp","")
if os.path.exists(optfolder + "/rating_table/ratTabIn.txt"):
os.remove(optfolder + "/rating_table/ratTabIn.txt")
if os.path.exists(optfolder + "/rating_table/userout.txt"):
os.remove(optfolder + "/rating_table/userout.txt")
if os.path.exists(optfolder + "/rating_table/rattabout.txt"):
os.remove(optfolder + "/rating_table/rattabout.txt")
#******************************************************************************************************
# open "rating table output" file in text editor
# [No. of Xsections = No. of Reaches]
#******************************************************************************************************
fn_lst=[]
tn_lst=[]
subriver = optfolder + "/subrivers.shp"
sr = arcpy.SearchCursor(subriver, "", "", "FROM_NODE;TO_NODE", "")
for node in sr:
fn = node.getValue("FROM_NODE")
tn = node.getValue("TO_NODE")
fn_lst.append(fn)
tn_lst.append(tn)
# use reach numbers to name "userout" and "rattabout" files
reach_lst = [x for x in fn_lst if x in tn_lst]
#******************************************************************************************************
# turn Xeditor OFF and Precipitation Depths ON
#******************************************************************************************************
# To make sure we have all the rating table text files for WinTR20 input text file
# Please make sure that text files pop up after transect tool have values NOT blank
myPath = optfolder + "/rating_table/"
txtCounter = len(glob.glob1(myPath,"*.txt"))
reachCounter = 2*len(reach_lst)
if txtCounter == reachCounter:
tool6.enabled = False
button13.enabled = True
else:
button13.enabled = False
self.Show(False) # close dialog box after deleting files and activating precip depth button
def OnRecalculate(self, event):
arcpy.env.scratchWorkspace = scratchfolder
arcpy.env.workspace = optfolder
#******************************************************************************************************
# get reach slope, Wbf, Dbf, Manning numbers, reach number, stream it drains to number,
# distance, and elevation to include in rating table input text file
#******************************************************************************************************
rs_val = self.reachslope.GetValue()
wbf_val = self.Wbf.GetValue()
dbf_val = self.Dbf.GetValue()
nMain_val = self.nMainValue.GetValue()
nLeft_val = self.nLeftValue.GetValue()
nRight_val = self.nRightValue.GetValue()
global rating
# use reach numbers to name "userout" and "rattabout" files
in_target = optfolder + "/trans_line.shp"
in_join = optfolder + "/subrivers.shp"
out_feature_class = optfolder + "/aux_folder/intersect_aux.shp"
arcpy.SpatialJoin_analysis(in_target,in_join,out_feature_class)
sr = arcpy.SearchCursor(out_feature_class, "", "", "FROM_NODE", "")
for node in sr:
rating = int(node.getValue("FROM_NODE"))
#******************************************************************************************************
# Rating table input text file string variables
#******************************************************************************************************
datastring = ""
datastring = datastring + "%s" "\n" %(rs_val)
datastring = datastring + "%s" "\n" %(wbf_val)
datastring = datastring + "%s" "\n" %(dbf_val)
datastring = datastring + "%s" "\n" %(nMain_val)
datastring = datastring + "%s" "\n" %(nLeft_val)
datastring = datastring + "%s" "\n" %(nRight_val)
datastring = datastring + "XS%s" "\n" %(rating)
datastring = datastring + "Reach%s" "\n" %(rating)
datastring = datastring + "\n"
x=""
for a,b in zip(trans_distance,trans_dem):
x = x + "{:<6}{:>7}".format(a,b)
x = x + "\n"
datastring = datastring + x
#******************************************************************************************************
# run executable using above input file to generate output file
# (input file have to be generated before running executable)
#******************************************************************************************************
if os.path.exists(optfolder + "/rating_table"):
path = optfolder + "/rating_table"
os.chdir(path)
defFN = optfolder + "/rating_table/ratTabIn.txt"
ratfile = open(defFN,"w")
ratfile.write(datastring)
ratfile.close()
os.system(path + "/rattab.exe")
time.sleep(4)
else:
path = optfolder + "/rating_table"
os.mkdir(path, 0755)
src = r"" + Directory + "/data/rating_table/rattab.exe"
shutil.copy2(src, path)
os.chdir(path)
defFN = optfolder + "/rating_table/ratTabIn.txt"
ratfile = open(defFN,"w")
ratfile.write(datastring)
ratfile.close()
os.system(path + "/rattab.exe")
time.sleep(4)
#******************************************************************************************************
# open "rating table output" file in text editor
# [No. of Xsections = No. of Reaches]
#******************************************************************************************************
arcpy.env.addOutputsToMap = True
rattabout_new = optfolder + "/rating_table/rattabout_reach" + str(rating) + ".txt"
rattabout_orig = optfolder + "/rating_table/rattabout.txt"
shutil.copy(rattabout_orig,rattabout_new)
userout_orig = optfolder + "/rating_table/userout.txt"
userout_new = optfolder + "/rating_table/userout_reach" + str(rating) + ".txt"
arcpy.CopyFeatures_management(optfolder + "/transect.shp", optfolder + "/xpoint_reach" + str(rating) + ".shp")
arcpy.CopyFeatures_management(optfolder + "/line_trun.shp", optfolder + "/xline_reach" + str(rating) + ".shp")
shutil.copy(userout_orig,userout_new)
os.system(optfolder + "/rating_table/userout_reach" + str(rating) + ".txt")
#******************************************************************************************************
# 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 == "xpoint_reach" + str(rating):
lyr.visible = False
self.SetWindowStyle( style | wx.STAY_ON_TOP )
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
user_prompt = pythonaddins.MessageBox("Do you want to continue or recalculate?", "Transect Line", 6)
if user_prompt == "Continue":
self.OnClose(event)
def OnPlot(self, event):
self.Bind(wx.EVT_CLOSE, self.OnClose) # 06-17-2015: Plot cross-section window close event
## get elevation and flow accumulation values along transect line
elevList = []
flowList = []
path = optfolder + "/transect.shp"
sc = arcpy.SearchCursor(path,"","","dem_x;flowacc_x","")
for att in sc:
elev = att.getValue("dem_x")
elevList.append(elev)
flow = att.getValue("flowacc_x")
flowList.append(flow)
wbf_val = self.Wbf.GetValue()
wbf = float(wbf_val)
dbf_val = self.Dbf.GetValue()
dbf = float(dbf_val)
demres = 49.2 # 1/29/2018: Horizontal distance is already in units of feet so this has to be in feet
plotlist = elevList
u = 0.65
# get minimum index value
minindx_elev = [idx for idx,val in enumerate(elevList) if val == min(elevList)]
minindx = int(minindx_elev[0]) # 04-03-2017 added to remove wiggle in the center
bfelev = min(elevList)
ileft = minindx
xcntr = demres * ileft + (wbf / 2)
alpha = (wbf / 2) / (dbf ** u)
numpts = int(len(elevList))
for i in range(minindx,numpts-2):
if elevList[i+1] < elevList[i]:
elevList[i+1] = elevList[i]
for j in range(minindx,1,-1):
if elevList[i-1] < elevList[i]:
elevList[i-1] = elevList[i]
ylist = []
xlist = []
for i in range(1,numpts+9):
ylist.append(0)
xlist.append(0)
for i in range(1,minindx+1):
ylist[i-1] = elevList[i-1]
xlist[i-1] = demres*(i-1)
for i in range(1,numpts-minindx):
ylist[(numpts+8)-i] = elevList[numpts-i]
xlist[(numpts+8)-i] = demres*(numpts-i)
ylist[minindx+4] = bfelev - dbf
xlist[minindx+4] = xcntr
for i in range(1,5):
indx1 = (minindx + 4) - i
indx2 = (minindx + 4) + i
ylist[indx1] = (bfelev - dbf) + ((dbf/4) * i)
ylist[indx2] = (bfelev - dbf) + ((dbf/4) * i)
xlist[indx1] = xcntr - (alpha * (((dbf/4) * i) ** u))
xlist[indx2] = xcntr + (alpha * (((dbf/4) * i) ** u))
numplotpoints = numpts + 8
ylmax = ylist[0]
yrmax = ylist[numplotpoints - 1]
ymax = min(ylmax,yrmax)
for i in range(minindx-1,1,-1):
if plotlist[i-1] >= ymax and plotlist[i] < ymax:
break
gplotList = ylist
minindx = [idx for idx,val in enumerate(gplotList) if val == min(gplotList)]
minindx = minindx[0]
iminmod = 0
for i in range((minindx - 1), 1, -1):
if ylist[i-1] >= ymax and ylist[i] < ymax:
iminmod = i
break
for i in range(minindx, numpts):
if plotlist[i-1] < ymax and plotlist[i] >= ymax:
break
for i in range(minindx - 1,numplotpoints):
if ylist[i-1] < ymax and ylist[i] >= ymax:
imaxmod = i + 1
break
xmod = []
ymod = []
for i in range(iminmod, imaxmod):
x = xlist[i-1]
xmod.append(x)
y = ylist[i-1]
ymod.append(y)
self.data = zip(xmod,ymod)
# "self.data" holds both x and y lists for plot
frm = wx.Frame(self, -1, "Elevation stage profile", size=(600,450))
client = plot.PlotCanvas(frm)
line = plot.PolyLine(self.data, legend="", colour="red", width=1)
gc = plot.PlotGraphics([line], "Height Cross-section", "Station (ft)", "Elevation (ft)")
client.Draw(gc, xAxis= (min(xmod)-20,max(xmod)+20), yAxis= (min(ymod)-10,max(ymod)+10))
frm.Show(True)
def OnExport(self, event):
elevList = []
flowList = []
path = optfolder + "/transect.shp"
sc = arcpy.SearchCursor(path,"","","dem_x;flowacc_x","")
for att in sc:
elev = att.getValue("dem_x")
elevList.append(elev)
flow = att.getValue("flowacc_x")
flowList.append(flow)
wbf_val = self.Wbf.GetValue()
wbf = float(wbf_val)
dbf_val = self.Dbf.GetValue()
dbf = float(dbf_val)
demres = 49.2 # 1/29/2018: Horizontal distance is already in units of feet so this has to be in feet
plotlist = elevList
u = 0.65
# get minimum index value
minindx_elev = [idx for idx,val in enumerate(elevList) if val == min(elevList)]
minindx = int(minindx_elev[0]) # 04-03-2017 added to remove wiggle in the center
bfelev = min(elevList)
ileft = minindx
xcntr = demres * ileft + (wbf / 2)
alpha = (wbf / 2) / (dbf ** u)
numpts = int(len(elevList))
for i in range(minindx,numpts-2):
if elevList[i+1] < elevList[i]:
elevList[i+1] = elevList[i]
for j in range(minindx,1,-1):
if elevList[i-1] < elevList[i]:
elevList[i-1] = elevList[i]
ylist = []
xlist = []
for i in range(1,numpts+9):
ylist.append(0)
xlist.append(0)
for i in range(1,minindx+1):
ylist[i-1] = elevList[i-1]
xlist[i-1] = demres*(i-1)
for i in range(1,numpts-minindx):
ylist[(numpts+8)-i] = elevList[numpts-i]
xlist[(numpts+8)-i] = demres*(numpts-i)
ylist[minindx+4] = bfelev - dbf
xlist[minindx+4] = xcntr
for i in range(1,5):
indx1 = (minindx + 4) - i
indx2 = (minindx + 4) + i
ylist[indx1] = (bfelev - dbf) + ((dbf/4) * i)
ylist[indx2] = (bfelev - dbf) + ((dbf/4) * i)
xlist[indx1] = xcntr - (alpha * (((dbf/4) * i) ** u))
xlist[indx2] = xcntr + (alpha * (((dbf/4) * i) ** u))
numplotpoints = numpts + 8
ylmax = ylist[0]
yrmax = ylist[numplotpoints - 1]
ymax = min(ylmax,yrmax)
for i in range(minindx-1,1,-1):
if plotlist[i-1] >= ymax and plotlist[i] < ymax:
break
gplotList = ylist
minindx = [idx for idx,val in enumerate(gplotList) if val == min(gplotList)]
minindx = minindx[0]
for i in range((minindx - 1), 1, -1):
if ylist[i-1] >= ymax and ylist[i] < ymax:
iminmod = i
break
for i in range(minindx, numpts):
if plotlist[i-1] < ymax and plotlist[i] >= ymax:
break
for i in range(minindx - 1,numplotpoints):
if ylist[i-1] < ymax and ylist[i] >= ymax:
imaxmod = i + 1
break
xmod = []
ymod = []
for i in range(iminmod, imaxmod):
x = xlist[i-1]
xmod.append(x)
y = ylist[i-1]
ymod.append(y)
xmod = ["%0.2f" %i for i in xmod]
ymod = ["%0.2f" %i for i in ymod]
# obtain reach number
reachpath = optfolder + "/reachslope.shp"
reachslope = arcpy.SearchCursor(reachpath, "", "", "FID_subshe", "")
for fid in reachslope:
val_fid1 = fid.getValue("FID_subshe")
subriverpath = optfolder + "/subrivers.shp"
subriver = arcpy.SearchCursor(subriverpath, "", "", "FID;GRID_CODE", "")
for gc in subriver:
val_fid2 = gc.getValue("FID")
if val_fid1 == val_fid2:
val_gc = gc.getValue("GRID_CODE")
#******************************************************************************************************
# open "rating table output" file in text editor
# [No. of Xsections = No. of Reaches]
#******************************************************************************************************
fn_lst=[]
tn_lst=[]
subriver = optfolder + "/subrivers.shp"
sr = arcpy.SearchCursor(subriver, "", "", "FROM_NODE;TO_NODE", "")
for node in sr:
fn = node.getValue("FROM_NODE")
tn = node.getValue("TO_NODE")
fn_lst.append(fn)
tn_lst.append(tn)
# use reach numbers to name "userout" and "rattabout" files
reach_lst = [x for x in fn_lst if x in tn_lst]
for r in zip(reach_lst):
if val_gc in r:
rating = val_gc
#******************************************************************************************************
# Write elevation stage profile and distance to text file
#******************************************************************************************************
datastring = ""
for a,b in zip(xmod,ymod):
datastring = datastring + "{:<6}{:>7}".format(a,b)
datastring = datastring + "\n"
if os.path.exists(optfolder + "/elev_stage_profile"):
elev_stage = optfolder + "/elev_stage_profile/elev_stage_reach" + str(rating) + ".txt"
elev_stage_file = open(elev_stage,"w")
elev_stage_file.write(datastring)
elev_stage_file.close()
else:
path = optfolder + "/elev_stage_profile"
os.mkdir(path, 0755)
elev_stage = optfolder + "/elev_stage_profile/elev_stage_reach" + str(rating) + ".txt"
elev_stage_file = open(elev_stage,"w")
elev_stage_file.write(datastring)
elev_stage_file.close()
pythonaddins.MessageBox("Elevation stage profile is successfully exported to %s" %elev_stage,"Elevation Stage Proile",0)
def OnApply(self, event):
arcpy.env.scratchWorkspace = scratchfolder
arcpy.env.workspace = optfolder
#******************************************************************************************************
# get reach slope, Wbf, Dbf, Manning numbers, reach number, stream it drains to number,
# distance, and elevation to include in rating table input text file
#******************************************************************************************************
rs_val = self.reachslope.GetValue()
wbf_val = self.Wbf.GetValue()
dbf_val = self.Dbf.GetValue()
nMain_val = self.nMainValue.GetValue()
nLeft_val = self.nLeftValue.GetValue()
nRight_val = self.nRightValue.GetValue()
global rating
in_target = optfolder + "/trans_line.shp"
in_join = optfolder + "/subrivers.shp"
out_feature_class = optfolder + "/aux_folder/intersect_aux.shp"
arcpy.SpatialJoin_analysis(in_target,in_join,out_feature_class)
sr = arcpy.SearchCursor(out_feature_class, "", "", "FROM_NODE", "")
for node in sr:
rating = int(node.getValue("FROM_NODE"))
#******************************************************************************************************
# Rating table input text file string variables
#******************************************************************************************************
datastring = ""
datastring = datastring + "%s" "\n" %(rs_val)
datastring = datastring + "%s" "\n" %(wbf_val)
datastring = datastring + "%s" "\n" %(dbf_val)
datastring = datastring + "%s" "\n" %(nMain_val)
datastring = datastring + "%s" "\n" %(nLeft_val)
datastring = datastring + "%s" "\n" %(nRight_val)
datastring = datastring + "XS%s" "\n" %(rating)
datastring = datastring + "Reach%s" "\n" %(rating)
datastring = datastring + "\n"
x=""
for a,b in zip(trans_distance,trans_dem):
x = x + "{:<6}{:>7}".format(a,b)
x = x + "\n"
datastring = datastring + x
#******************************************************************************************************
# run executable using above input file to generate output file
# (input file have to be generated before running executable)
#******************************************************************************************************
if os.path.exists(optfolder + "/rating_table"):
path = optfolder + "/rating_table"
os.chdir(path)
defFN = optfolder + "/rating_table/ratTabIn.txt"
ratfile = open(defFN,"w")
ratfile.write(datastring)
ratfile.close()
os.system(path + "/rattab.exe")
time.sleep(4)
else:
path = optfolder + "/rating_table"
os.mkdir(path, 0755)
src = r"" + Directory + "/data/rating_table/rattab.exe"
shutil.copy2(src, path)
os.chdir(path)
defFN = optfolder + "/rating_table/ratTabIn.txt"
ratfile = open(defFN,"w")
ratfile.write(datastring)
ratfile.close()
os.system(path + "/rattab.exe")
time.sleep(4)
#******************************************************************************************************
# open "rating table output" file in text editor
# [No. of Xsections = No. of Reaches]
#******************************************************************************************************
rattabout_new = optfolder + "/rating_table/rattabout_reach" + str(rating) + ".txt"
rattabout_orig = optfolder + "/rating_table/rattabout.txt"
shutil.copy(rattabout_orig,rattabout_new)
userout_orig = optfolder + "/rating_table/userout.txt"
userout_new = optfolder + "/rating_table/userout_reach" + str(rating) + ".txt"
arcpy.CopyFeatures_management(optfolder + "/transect.shp", optfolder + "/xpoint_reach" + str(rating) + ".shp")
arcpy.CopyFeatures_management(optfolder + "/line_trun.shp", optfolder + "/xline_reach" + str(rating) + ".shp")
shutil.copy(userout_orig,userout_new)
os.system(optfolder + "/rating_table/userout_reach" + str(rating) + ".txt")
#******************************************************************************************************
# 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 == "xpoint_reach" + str(rating):
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "intersect_aux":
arcpy.mapping.RemoveLayer(df, lyr)
if lyr.name == "trans_line":
arcpy.mapping.RemoveLayer(df, lyr)
arcpy.Delete_management(optfolder + "/aux_folder/intersect_aux.shp")
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
self.btnApply.Disable() # 1/3/2017: Disable "Apply" button after first execution
self.btnRecalculate.Enable() # 1/3/2017: activate "Recalculate" option once first run of rating table is executed
user_prompt = pythonaddins.MessageBox("Do you want to continue or recalculate?", "Transect Line", 6)
if user_prompt == "Continue":
self.OnClose(event)