GISHydroNXT System Documentation

#*********************************************************************************
# Author:       UMD
# Date:         24-07-2018
# Modified:     n/a
# Classes:      ThomasDischarge()
# Functions:    hydro.CheckGages() ; hydro.FC() ; hydro.pctAsoil() ; hydro.pctCsoil()
                hydro.pctDsoil(); hydro.FloodIntervalLists() ; hydro.TaskerString()
                hydro.openbrowser()
# Modules:      arcpy ; os ; webbrowser
# Comments:     n/a
#*********************************************************************************
class ThomasDischarge(object):
    """Implementation for GISHydroNXT_addin.button4 (Button)"""
    def __init__(self):
        self.enabled = False
        self.checked = False
    def onClick(self):
        arcpy.env.scratchWorkspace = scratchfolder
        arcpy.env.workspace = optfolder
        #*******************************************************************************************************
        # gage checking and selection
        #*******************************************************************************************************
        # added on 10-13-2017: for to include in intersect analysis
        usgsgages = r"" + Directory + "/data/maryland/usgsgagesm.shp"
        mdgages = r"" + Directory + "/data/maryland/mdgagedstreams2016.shp"
        outletpoint = optfolder + "/outletpoint.shp"
        mask = optfolder + "/mask.shp"
        gagefound = hydro.CheckGages(usgsgages,mdgages,outletpoint,mask,optfolder)

        global gagelist
        gagelist = gagefound

        #*******************************************************************************************************
        # Check if gauge exist in the "gagefound" list and implement if/else for
        # gage dialog to pop up
        #*******************************************************************************************************
        if not gagelist:
            #*****************************************************************************
            # read global variables (defined in basin stat) into variable names for use
            # in Thomas peak discharge analysis
            FC = float(hydro.FC) # LI is already declared as global variable
            FC = "{0:.2f}".format(FC)
            DA = areami2
            HA = float(hydro.pctAsoil)
            HC = float(hydro.pctCsoil)
            HD = float(hydro.pctDsoil)
            SLL = float(landslope)
            HCD = float(HC+HD)
            ImpA = float(IA)

            #*****************************************************************************
            # Define initial values and lists. add and index those lists as part of
            # dictionary
            #*****************************************************************************
            sumarea = 0
            provstring = ""
            outtaskerstring = ""
            sQ1p25 = 0
            sQ1p50 = 0
            sQ1p75 = 0
            sQ2 = 0
            sQ5 = 0
            sQ10 = 0
            sQ25 = 0
            sQ50 = 0
            sQ100 = 0
            sQ200 = 0
            sQ500 = 0
            Q1p25list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q1p50list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q2list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q5list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q10list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q25list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q50list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q100list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q200list = [0, 0, 0, 0, 0, 0, 0, 0]
            Q500list = [0, 0, 0, 0, 0, 0, 0, 0]

            qlist = {"1":[],"2":[],"3":[],"4":[],"5":[],"6":[],"7":[],"8":[],"9":[],"10":[]}
            qlist["1"].extend(Q1p25list)
            qlist["2"].extend(Q1p50list)
            qlist["3"].extend(Q2list)
            qlist["4"].extend(Q5list)
            qlist["5"].extend(Q10list)
            qlist["6"].extend(Q25list)
            qlist["7"].extend(Q50list)
            qlist["8"].extend(Q100list)
            qlist["9"].extend(Q200list)
            qlist["10"].extend(Q500list)

            #*****************************************************************************
            # copy root tasker files into optfolder 
            #*****************************************************************************
            src = r"" + Directory + "/data/tasker"
            dst = optfolder + "/tasker"
            shutil.copytree(src, dst)

            #*****************************************************************************
            # read zonal stat table (theVTab), declare fields into variables, and
            # loop through theVTab fields to getValue
            #*****************************************************************************
            theVTab = optfolder + "/theVTab.dbf"
            theVTab = arcpy.SearchCursor("theVTab","","","Count","")

            #*****************************************************************************
            # loop to get total count of pixels of basingrid
            #*****************************************************************************
            for each in theVTab:
                count = each.getValue("Count")
                sumarea = sumarea+count
            sumArea = sumarea

            del each

            theVTab = optfolder + "/theVTab.dbf"
            theVTab = arcpy.SearchCursor("theVTab","","","Province;Count;Q1.25;Q1.50;Q1.75;Q2;Q5;Q10;Q25;Q50;Q100;Q200;Q500","")

            # loop to begin tasker handling and area weighted analysis
            discharge_values = []
            flood_intervals = []
            for row in theVTab:
                AreaField = float(row.getValue("Count"))
                areapercent = float((AreaField/sumArea)* 100)
                if row.getValue("Province") == "A":
                    provstring = provstring + "       -Appalachian Plateaus and Allegheny Ridges %s percent of area" "\n" %("{0:.2f}".format(areapercent))
                elif row.getValue("Province") == "B":
                    provstring = provstring + "       -Blue Ridge and Great Valley %s percent of area" "\n" %("{0:.2f}".format(areapercent))
                elif row.getValue("Province") == "P":
                    provstring = provstring + "       -Piedmont %s percent of area" "\n" %("{0:.2f}".format(areapercent))
                elif row.getValue("Province") == "W":
                    provstring = provstring + "       -Western Coastal Plain %s percent of area" "\n" %("{0:.2f}".format(areapercent))
                elif row.getValue("Province") == "E":
                    provstring = provstring + "       -Eastern Coastal Plain %s percent of area" "\n" %("{0:.2f}".format(areapercent))

                else:
                    pythonaddins.MessageBox("No Province Selected", "Problem...", 0)

                intaskerstring = "thomasout.txt" "\n"
                intaskerstring = intaskerstring + "" "\n"
                if row.getValue("Province") == "A":
                    intaskerstring = intaskerstring + "a"  "\n"
                    intaskerstring = intaskerstring + "%s" "\n" %(DA)
                    intaskerstring = intaskerstring + "%s" "\n" %(SLL)
                elif row.getValue("Province") == "B":
                    intaskerstring = intaskerstring + "p"  "\n"
                    intaskerstring = intaskerstring + "%s" "\n" %(DA)
                    intaskerstring = intaskerstring + "%s" "\n" %(FC)
                    intaskerstring = intaskerstring + "%s" "\n" %(LI)
                    intaskerstring = intaskerstring + "%s" "\n" %(ImpA)
                elif row.getValue("Province") == "P":
                    intaskerstring = intaskerstring + "p"  "\n"
                    intaskerstring = intaskerstring + "%s" "\n" %(DA)
                    intaskerstring = intaskerstring + "%s" "\n" %(FC)
                    intaskerstring = intaskerstring + "%s" "\n" %(LI)
                    intaskerstring = intaskerstring + "%s" "\n" %(ImpA)

                elif row.getValue("Province") == "W":
                    intaskerstring = intaskerstring + "wc" "\n"
                    intaskerstring = intaskerstring + "%s" "\n" %(DA) 
                    intaskerstring = intaskerstring + "%s" "\n" %(ImpA)
                    intaskerstring = intaskerstring + "%s" "\n" %(HCD)
                elif row.getValue("Province") == "E":
                    intaskerstring = intaskerstring + "ec" "\n"
                    intaskerstring = intaskerstring + "%s" "\n" %(DA)
                    intaskerstring = intaskerstring + "%s" "\n" %(SLL)
                    intaskerstring = intaskerstring + "%s" "\n" %(HA)

                intaskerstring = intaskerstring + "N" "\n"
                intaskerstring = intaskerstring + "N" "\n"
                thomas2016 = optfolder + "/tasker/thomas2016.txt"
                thomasin = open(thomas2016, "w")
                thomasin.write(intaskerstring)
                thomasin.close()

                os.chdir(optfolder + "/tasker")
                os.system(optfolder + "/tasker/thomas2016auto.exe")
                time.sleep(4)

                #*****************************************************************************
                # open "thomasout.txt" file and read peak discharges
                #*****************************************************************************
                infilename = optfolder + "/tasker/thomasout.txt"
                infile = open(infilename,"r").readlines()

                for i in infile:
                    outtaskerstring = outtaskerstring + i

                #*****************************************************************************
                # Discharge values for each province for "theVTab" are not possible to write
                # at this time because update cursor can"t be used within search cursor "for"
                # loop -- appended to list "discharge_values" and will be used at the end to
                # update discharge values in "theVTab" 
                #*****************************************************************************
                theline = infile[11].split()
                Q1p25 = theline[1]
                discharge_values.append(Q1p25)
                theline = infile[12].split()
                Q1p50 = theline[1]
                discharge_values.append(Q1p50)
                Q1p75 = -999
                theline = infile[13].split()
                Q2 = theline[1]
                discharge_values.append(Q2)
                theline = infile[14].split()
                Q5 = theline[1]
                discharge_values.append(Q5)
                theline = infile[15].split()
                Q10 = theline[1]
                discharge_values.append(Q10)
                theline = infile[16].split()
                Q25 = theline[1]
                discharge_values.append(Q25)
                theline = infile[17].split()
                Q50 = theline[1]
                discharge_values.append(Q50)
                theline = infile[18].split()
                Q100 = theline[1]
                discharge_values.append(Q100)
                theline = infile[19].split()
                Q200 = theline[1]
                discharge_values.append(Q200)
                theline = infile[20].split()
                Q500 = theline[1]
                discharge_values.append(Q500)

                #*****************************************************************************
                # compute discharge and assign confidence intervals to qlist entries
                #*****************************************************************************
                Q1p25 = float(Q1p25)
                sQ1p25 = sQ1p25 + (Q1p25* AreaField)
                Q1p50 = float(Q1p50)
                sQ1p50 = sQ1p50 + (Q1p50* AreaField)
                Q1p75 = float(Q1p75)
                sQ1p75 = sQ1p75 + (Q1p75* AreaField)
                Q2 = float(Q2)
                sQ2 = sQ2 + (Q2* AreaField)
                Q5 = float(Q5)
                sQ5 = sQ5 + (Q5* AreaField)
                Q10 = float(Q10)
                sQ10 = sQ10 + (Q10* AreaField)
                Q25 = float(Q25)
                sQ25 = sQ25 + (Q25* AreaField)
                Q50 = float(Q50)
                sQ50 = sQ50 + (Q50* AreaField)
                Q100 = float(Q100)
                sQ100 = sQ100 + (Q100* AreaField)
                Q200 = float(Q200)
                sQ200 = sQ200 + (Q200* AreaField)
                Q500 = float(Q500)
                sQ500 = sQ500 + (Q500* AreaField)

                for i, lines in enumerate(infile):
                    if i > 24 and i < 35:
                        line = lines.split()
                        a1 = line[1]
                        a2 = line[2]
                        a3 = line[3]
                        a4 = line[4]
                        a5 = line[5]
                        a6 = line[6]
                        a7 = line[7]
                        a8 = line[8]
                        c1 = [(float(a1)* areapercent)/100]
                        c2 = [(float(a2)* areapercent)/100]
                        c3 = [(float(a3)* areapercent)/100]
                        c4 = [(float(a4)* areapercent)/100]
                        c5 = [(float(a5)* areapercent)/100]
                        c6 = [(float(a6)* areapercent)/100]
                        c7 = [(float(a7)* areapercent)/100]
                        c8 = [(float(a8)* areapercent)/100]
                        qlist = [c1,c2,c3,c4,c5,c6,c7,c8]
                        flood_intervals.append(qlist)

            
            lists = {i:[el[0] for el in v] for i, v in enumerate(flood_intervals, start=1)}
                    
            del row

            #*****************************************************************************
            # index and count number of sub-lists -- prepare for TaskerString function
            #*****************************************************************************
            number = hydro.FloodIntervalLists(lists)
            if number > 10:
                q_list1 = [sum(i) for i in zip(lists[1],lists[11])]
                q_list2 = [sum(i) for i in zip(lists[2],lists[12])]
                q_list3 = [sum(i) for i in zip(lists[3],lists[13])]
                q_list4 = [sum(i) for i in zip(lists[4],lists[14])]
                q_list5 = [sum(i) for i in zip(lists[5],lists[15])]
                q_list6 = [sum(i) for i in zip(lists[6],lists[16])]
                q_list7 = [sum(i) for i in zip(lists[7],lists[17])]
                q_list8 = [sum(i) for i in zip(lists[8],lists[18])]
                q_list9 = [sum(i) for i in zip(lists[9],lists[19])]
                q_list10 = [sum(i) for i in zip(lists[10],lists[20])]
            else:
                q_list1 = lists[1]
                q_list2 = lists[2]
                q_list3 = lists[3]
                q_list4 = lists[4]
                q_list5 = lists[5]
                q_list6 = lists[6]
                q_list7 = lists[7]
                q_list8 = lists[8]
                q_list9 = lists[9]
                q_list10 = lists[10]

            #*****************************************************************************
            # discharge computation based on province
            #*****************************************************************************
            Q1p25 = int(sQ1p25/sumArea)
            Q1p50 = int(sQ1p50/sumArea)
            Q1p75 = int(sQ1p75/sumArea)
            Q2 = int(sQ2/sumArea)
            Q5 = int(sQ5/sumArea)
            Q10 = int(sQ10/sumArea)
            Q25 = int(sQ25/sumArea)
            Q50 = int(sQ50/sumArea)
            Q100 = int(sQ100/sumArea)
            Q200 = int(sQ200/sumArea)
            Q500 = int(sQ500/sumArea)

            #*****************************************************************************
            # Basin statistics analysis date
            #*****************************************************************************
            now = datetime.now()
            month = now.strftime("%B")
            day = now.strftime("%d")
            year = now.strftime("%Y")

            #*****************************************************************************
            # Text file string variables
            #*****************************************************************************
            datastring = ""
            datastring = datastring + "GISHydro Release Version Date: %s" "\n" %(Modifieddt)
            datastring = datastring + "Project Name:                  %s" %(proj)
            datastring = datastring + "" "\n"
            datastring = datastring + "Analysis Date:                 %s %s, %s " "\n" %(month,day,year) 
            datastring = datastring + "Geographic Province(s):" "\n"
            datastring = datastring + provstring
            datastring = datastring + "" "\n"
            datastring = datastring + "Q(1.25):   %s cfs" "\n" %(Q1p25)
            datastring = datastring + "Q(1.50):   %s cfs" "\n" %(Q1p50)
            datastring = datastring + "Q(2):      %s cfs" "\n" %(Q2)
            datastring = datastring + "Q(5):      %s cfs" "\n" %(Q5)
            datastring = datastring + "Q(10):     %s cfs" "\n" %(Q10)
            datastring = datastring + "Q(25):     %s cfs" "\n" %(Q25)
            datastring = datastring + "Q(50):     %s cfs" "\n" %(Q50)
            datastring = datastring + "Q(100):    %s cfs" "\n" %(Q100)
            datastring = datastring + "Q(200):    %s cfs" "\n" %(Q200)
            datastring = datastring + "Q(500):    %s cfs" "\n" %(Q500)

            datastring = datastring + "" "\n"
            datastring = datastring + "Area Weighted Prediction Intervals (from Tasker)" "\n"
            datastring = datastring + " Return     50 PERCENT        67 PERCENT        90 PERCENT        95 PERCENT" "\n"
            datastring = datastring + " Period  lower    upper    lower    upper    lower    upper    lower    upper" "\n"
            datastring = datastring + hydro.TaskerString(1.25,q_list1)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(1.50,q_list2)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(2,q_list3)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(5,q_list4)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(10,q_list5)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(25,q_list6)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(50,q_list7)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(100,q_list8)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(200,q_list9)
            datastring = datastring + "" "\n"
            datastring = datastring + hydro.TaskerString(500,q_list10)
            datastring = datastring + "" "\n"

            datastring = datastring + "" "\n"
            datastring = datastring + "" "\n"
            datastring = datastring + "Individual Province Tasker Analyses Follow: " "\n"
            datastring = datastring + ""
            datastring = datastring + outtaskerstring

            #*****************************************************************************
            # write strings to basin stat text file.
            #*****************************************************************************
            defFN = optfolder + "/frdischarges.txt"
            fr = open(defFN,"w")
            fr.write(datastring)
            fr.close()

            #*****************************************************************************
            # open "frdischarges" file in text editor
            #*****************************************************************************
            hydro.openbrowser(defFN)

            #*****************************************************************************
            # 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 == "mdlayer":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "gagefield":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoint":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoly":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "limepoly":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoint":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletbuffer":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "gagefield":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoly":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "mask_ints": # added on 10-22-2017: this layer is output of intersect tool
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "gauge_outlet": # added on 10-22-2017: this layer is output of join and intersect tool
                    arcpy.mapping.RemoveLayer(df, lyr)

            arcpy.RefreshTOC()
            arcpy.RefreshActiveView()

            #*****************************************************************************
            # turn peak discharge OFF, and "S" & Add Streams ON
            #*****************************************************************************
            tool7.enabled = False
            button4.enabled = False
            button5_1.enabled = True
            button6.enabled = True
            tool4.enabled = True
        else:
            mxd = arcpy.mapping.MapDocument("CURRENT")
            df = arcpy.mapping.ListDataFrames(mxd)[0]
            layers = arcpy.mapping.ListLayers(mxd, "", df)
            for lyr in layers:
                if lyr.name == "mdlayer":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "gagefield":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoint":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoly":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "limepoly":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoint":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletbuffer":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "gagefield":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "outletpoly":
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "mask_ints": # added on 10-22-2017: this layer is output of intersect tool
                    arcpy.mapping.RemoveLayer(df, lyr)
                if lyr.name == "gauge_outlet": # added on 10-22-2017: this layer is output of join and intersect tool
                    arcpy.mapping.RemoveLayer(df, lyr)
            dlg2 = GageListThomas()