GISHydroNXT System Documentation

#*********************************************************************************
# Author:       UMD
# Date:         24-07-2018
# Modified:     n/a
# Classes:      DelineateSubwatersheds() 
# Functions:    n/a
# Modules:      arcpy ; os
# Comments:     n/a
#*********************************************************************************
class DelineateSubwatersheds(object):
    """Implementation for GISHydroNXT_addin.button8 (Button)"""
    def __init__(self):
        self.enabled = False
        self.checked = False
    def onClick(self):
        # extent is smaller when subwatersheds are delineated again after Reset
        arcpy.env.extent = "MAXOF"
        arcpy.env.scratchWorkspace = scratchfolder
        arcpy.env.workspace = optfolder
        #*******************************************************************************************************        
        # add subwatersheds to view -- addOutputsToMap 
        #*******************************************************************************************************
        arcpy.env.addOutputsToMap = False
        basingrid = optfolder + "/basingrid"
        flowacc = optfolder + "/flowacc"
        flwdir = arcpy.sa.Times(optfolder + "/flowdir",basingrid)
        strlnk = arcpy.sa.StreamLink(optfolder + "/ModStreams",flwdir)
        strlnk.save(optfolder + "/strlnk")
        zonemax = arcpy.sa.ZonalStatistics(strlnk, "Value", flowacc, "MAXIMUM", "NODATA")
        faccmax = arcpy.sa.Con(flowacc == zonemax, zonemax, IsNull(zonemax))
        outlet = arcpy.sa.Con(faccmax > 0, faccmax)
        outlet.save(optfolder + "/outlets")

        # if outlets are added by user then add mod stream and user specified outlets 
        outlets_user   = optfolder + "/outlets_user"
        if os.path.exists(outlets_user):
            outlets_modStr = arcpy.sa.Con(arcpy.sa.IsNull(optfolder + "/outlets"), 0, optfolder + "/outlets")
            outlets_modStr.save(optfolder + "/outlets_str")
            arcpy.env.extent = optfolder + "/flowacc"
            arcpy.env.snapRaster = optfolder + "/flowacc"
            aux_outlet = arcpy.sa.SetNull(optfolder + "/outlets_str",1,"Value = 0")
            aux_outlet2 = arcpy.RasterToPoint_conversion(aux_outlet)
            aux_outlet3 = arcpy.Merge_management([optfolder + "/AddasOutlets.shp",aux_outlet2])
            subwshed = arcpy.sa.Watershed(optfolder + "/flowdir", aux_outlet3,"FID")
            subwshed.save(optfolder + "/subshed_temp")
            arcpy.env.extent = "MAXOF"
        else:
            subwshed = arcpy.sa.Watershed(optfolder + "/flowdir", optfolder + "/outlets", "VALUE")
            subwshed.save(optfolder + "/subshed_temp")

        arcpy.env.addOutputsToMap = True
        arcpy.RasterToPolygon_conversion(optfolder + "/subshed_temp",optfolder + "/tmpsubwshd","NO_SIMPLIFY","VALUE")
        tmpsub = optfolder + "/tmpsubwshd.shp"
        subshed = optfolder + "/subshed.shp"
        arcpy.Dissolve_management(tmpsub,subshed,"GRIDCODE","#","MULTI_PART","DISSOLVE_LINES")
        subrivers = optfolder + "/subrivers.shp"

        # new stream link raster to handle added outlets
        if os.path.exists(outlets_user):
            StrmFDRgrid = arcpy.sa.Times(optfolder + "/ModStreams",flwdir)
            newLnkGrid = arcpy.sa.Watershed(StrmFDRgrid, aux_outlet3, "FID")
            newLnkGrid.save(optfolder + "/newStrlnk")
            arcpy.sa.StreamToFeature(optfolder + "/newStrlnk", flwdir, subrivers, "NO_SIMPLIFY") # 06-23-2014: Modified to use new stream link raster"
            del aux_outlet2,aux_outlet3
        else:
            arcpy.sa.StreamToFeature(optfolder + "/ModStreams", flwdir, subrivers, "NO_SIMPLIFY")

        #*******************************************************************************************************
        # correct sub-basin indexing -- order of sub-basins has to be the same as
        # appearing in "subriver.shp"
        #*******************************************************************************************************
        target1 = optfolder + "/subrivers.shp"
        joint1  = optfolder + "/subshed.shp"
        spatial_join1 = optfolder + "/subriver_subshed.shp"

        # add table using spatial join (target = subriver, join = subshed)
        fieldmappings = arcpy.FieldMappings()
        fieldmappings.addTable(target1)
        fieldmappings.addTable(joint1)
        arcpy.SpatialJoin_analysis(target1, joint1, spatial_join1, "#", "#", fieldmappings, "HAVE_THEIR_CENTER_IN", "#", "#")

        # get subriver FID and centroid (or midpoint) XY values 
        sub_cur = arcpy.SearchCursor(spatial_join1, "", "", "Shape;ARCID;GRIDCODE", "")
        x = []
        y = []
        for row in sub_cur:
            XMidPoint = row.shape.positionAlongLine(0.50,True).firstPoint.X
            x.append(XMidPoint)
            YMidPoint = row.shape.positionAlongLine(0.50,True).firstPoint.Y
            y.append(YMidPoint)

        del sub_cur

        xy = list(zip(x,y))

        # create a point shapefile, add fields "ARCID" & "GRIDCODE", and add midpoint XY values to field "SHAPE@" 
        spatial_reference = arcpy.Describe(optfolder + "/subriver_subshed.shp").spatialReference
        arcpy.CreateFeatureclass_management(optfolder, "subriver_xy.shp", "POINT", "", "ENABLED", "DISABLED", spatial_reference)
        subriver_xy = optfolder + "/subriver_xy.shp"

        # added on 09-30-2014 to debug field type error -- changed field type to "Double" with precision "10"
        arcpy.AddField_management(subriver_xy,"ARCID","DOUBLE",10,"")
        arcpy.AddField_management(subriver_xy,"GRIDCODE","DOUBLE",10,"")
        sr_xy = arcpy.da.InsertCursor(subriver_xy, ("SHAPE@XY"))
        for i in xy:
            sr_xy.insertRow([i])

        # update "ARCID" and "GRIDCODE" row values using "subriver_subshed" shapefile
        sub_cur = arcpy.SearchCursor(spatial_join1, "", "", "ARCID;GRIDCODE", "")
        ic = arcpy.UpdateCursor(subriver_xy, "", "", "ARCID;GRIDCODE", "")

        row1 = sub_cur.next()
        row2 = ic.next()

        while row1:
            row2.setValue("ARCID", row1.getValue("ARCID"))
            ic.updateRow(row2)
            row2.setValue("GRIDCODE", row1.getValue("GRIDCODE"))
            ic.updateRow(row2)
            row1 = sub_cur.next()
            row2 = ic.next()

        del sub_cur, ic

        # perform spatial join (target = subshed, join = subriver_xy) and lists extraction
        target2 = optfolder + "/subshed.shp"
        joint2  = optfolder + "/subriver_xy.shp"
        spatial_join2  = optfolder + "/subshed_subriver.shp"

        fieldmappings = arcpy.FieldMappings()
        fieldmappings.addTable(target2)
        fieldmappings.addTable(joint2)
        arcpy.SpatialJoin_analysis(target2, joint2, spatial_join2, "#", "#", fieldmappings, "INTERSECT", "#", "#")

        # create a polygon which will contain subriver FID, ARCID, and GRIDCODE (obtained from spatially joined shapefile)
        spatial_reference = arcpy.Describe(optfolder + "/subshed_subriver.shp").spatialReference
        arcpy.CreateFeatureclass_management(optfolder, "subshed.shp", "POLYGON", "", "DISABLED", "DISABLED", spatial_reference)

        sc = arcpy.da.SearchCursor(spatial_join2, ("ARCID", "GRIDCODE", "SHAPE@"))

        storage = []
        for row in sc: # changed sF to sc [could be a mistake to have "sF"]
            storage.append(row)
            
        sortedlist = sorted(storage, key=lambda x:x[0]) # list sort is based on "ARCID"

        poly = optfolder + "/subshed.shp"

        # added on 09-30-2014 to debug field type error -- changed field type to "Double" with precision "10"
        arcpy.AddField_management(poly,"ARCID","DOUBLE",10,"")
        arcpy.AddField_management(poly,"GRIDCODE","DOUBLE",10,"")
        
        arcpy.DeleteField_management(poly,"Id")
        fcout = arcpy.da.InsertCursor(poly, ("ARCID","GRIDCODE","SHAPE@"))

        for i in sortedlist:
            fcout.insertRow(i)

        del fcout, sc

        # convert newly created and indexed sub-watershed shapefile to raster
        arcpy.PolygonToRaster_conversion(subshed, "ARCID", optfolder + "/subsheds", "", "ARCID", 30) # priority field is "ARCID"
        
        #*******************************************************************************************************
        # 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 == "tmpsubwshd":
                arcpy.mapping.RemoveLayer(df, lyr)
            if lyr.name == "subsheds":
                df.extent = lyr.getSelectedExtent()
                arcpy.mapping.RemoveLayer(df, lyr)
            if lyr.name == "subshed_subriver":
                arcpy.mapping.RemoveLayer(df, lyr)
            if lyr.name == "subriver_subshed":
                arcpy.mapping.RemoveLayer(df, lyr)
            if lyr.name == "subriver_xy":
                arcpy.mapping.RemoveLayer(df, lyr)
            if lyr.name == "subshed":
                arcpy.ApplySymbologyFromLayer_management(lyr,r"" + Directory + "/data/mdfiles/legends/watershed.lyr")
            if lyr.name == "subrivers":
                arcpy.ApplySymbologyFromLayer_management(lyr,r"" + Directory + "/data/mdfiles/legends/subrivers.lyr")
        arcpy.env.extent = "MAXOF"

        #*******************************************************************************************************
        # turn subwatershed delineation OFF
        #*******************************************************************************************************
        
        button7.enabled = False
        button8.enabled = False
        button9.enabled = True
        tool5.enabled = False