Implemented the Drainage area check as a check box

RCAT.pyt

@@ -257,8 +257,15 @@ class VBETtool(object):
             direction="Input")
         param20.value = 50000
 
+        param21 = arcpy.Parameter(
+            displayName="Validate Drainage Area Using ReachDist",
+            name="check_drain_area",
+            datatype="GPBoolean",
+            parameterType="Optional",
+            direction="Input")
+
         return [param0, param1, param2, param3, param4, param5, param6, param7, param8, param9, param10, param11,
-                param12, param13, param14, param15, param16, param17, param18, param19, param20]
+                param12, param13, param14, param15, param16, param17, param18, param19, param20, param21]
 
     def isLicensed(self):
         """Set whether tool is licensed to execute."""
@@ -298,7 +305,8 @@ class VBETtool(object):
                   p[17].valueAsText,
                   p[18].valueAsText,
                   p[19].valueAsText,
-                  p[20].valueAsText)
+                  p[20].valueAsText,
+                  p[21].valueAsText)
         return
 
 

RCAT_Drainage_Area_Check.py

@@ -0,0 +1,169 @@
+# -------------------------------------------------------------------------------
+# Name:        Drainage_Area_Check
+# Purpose:     Looks through the stream network for reaches that have lower drainage networks than reaches up stream of
+#              them, and modifies the network to fix that
+#
+# Author:      Braden Anderson
+#
+# Created:     06/2018
+# -------------------------------------------------------------------------------
+
+import arcpy
+import os
+from RCAT_Stream_Objects import DAValueCheckStream, ProblemStream, StreamHeap
+
+
+def main(stream_network):
+    """
+    The main function
+    :param stream_network: The stream network that we want to fix up
+    :return:
+    """
+    stream_heaps = find_streams(stream_network)
+    #check_heap(stream_network, stream_heaps)
+
+    problem_streams = find_problem_streams(stream_heaps)
+    #check_problem_streams(stream_network, problem_streams)
+
+    fix_problem_streams(stream_network, problem_streams)
+
+
+def find_streams(stream_network):
+    """
+    Creates a list of heaps, sorted by distance from the stream head
+    :param stream_network: The stream network to be used
+    :return:
+    """
+    arcpy.AddMessage("Finding Streams...")
+    stream_heaps = []
+    req_fields = ["ReachID", "StreamID", "ReachDist", "DA_sqkm"]
+    with arcpy.da.SearchCursor(stream_network, req_fields) as cursor:
+        for reach_id, stream_id, downstream_dist, drainage_area in cursor:
+            new_stream = DAValueCheckStream(reach_id, stream_id, downstream_dist, drainage_area)
+            new_stream_heap_index = find_new_stream_heap_index(stream_id, stream_heaps)
+
+            if new_stream_heap_index is not None:
+                stream_heaps[new_stream_heap_index].push_stream(new_stream)
+            else:
+                new_stream_heap = StreamHeap(new_stream)
+                stream_heaps.append(new_stream_heap)
+    return stream_heaps
+
+
+def check_heap(stream_network, stream_heaps):
+    with open(os.path.join(os.path.dirname(stream_network), "streams.txt"), 'w') as file:
+        for stream_heap in stream_heaps:
+            file.write(str(stream_heap) + '\n')
+    for stream_heap in stream_heaps:
+        streams = stream_heap.streams
+        for k in range(len(streams)):
+            try:
+                high_dist = streams[k].downstream_dist
+                low_dist_one = streams[(k*2) + 1].downstream_dist
+                low_dist_two = streams[(k*2) + 2].downstream_dist
+                if high_dist < low_dist_one or high_dist < low_dist_two:
+                    raise Exception("Error in stream id: " + str(streams[k].stream_id))
+            except IndexError:
+                pass
+    arcpy.AddMessage("Stream Heaps passed check")
+
+
+
+def find_new_stream_heap_index(stream_id, stream_heaps):
+    """
+    Finds the index of the heap that the stream belongs to
+    :param stream_id: The stream_id that we want to find in the heaps
+    :param stream_heaps: A list of heaps
+    :return: A number if that stream ID is in the list of heaps, otherwise, None
+    """
+    for i in range(len(stream_heaps)):
+        if stream_id == stream_heaps[i].stream_id:
+            return i
+    return None
+
+
+
+def find_problem_streams(stream_heaps):
+    """
+    Looks through the stream heaps, identifies streams that need to be fixed, and puts them in a list
+    :param stream_heaps: A list of stream heaps
+    :return:
+    """
+    arcpy.AddMessage("Identifying problem streams...")
+    problem_streams = []
+    for stream_heap in stream_heaps:
+        while len(stream_heap.streams) > 0:
+            downstream_reach = stream_heap.pop()
+            max_upstream_drainage_area = 0.0
+            for stream in stream_heap.streams:
+                if stream.drainage_area > max_upstream_drainage_area:
+                    max_upstream_drainage_area = stream.drainage_area
+
+            if downstream_reach.drainage_area < max_upstream_drainage_area:
+                new_problem_stream = ProblemStream(downstream_reach.reach_id, downstream_reach.stream_id, downstream_reach.drainage_area, max_upstream_drainage_area)
+                problem_streams.append(new_problem_stream)
+
+    return problem_streams
+
+
+def check_problem_streams(stream_network, problem_streams):
+    """
+    A simple function that's mean to write the data to a text file and raise errors if something seems wrong
+    """
+    max_orig_DA = 0
+    max_orig_DA_id = -1
+    with open(os.path.join(os.path.dirname(stream_network), "problemStreams.txt"), 'w') as file:
+        for problem_stream in problem_streams:
+            file.write(str(problem_stream) + '\n')
+            if problem_stream.orig_drainage_area > 50:
+                arcpy.AddWarning("Reach " + str(problem_stream.reach_id) + " may not be a problem stream")
+            if problem_stream.orig_drainage_area > max_orig_DA:
+                max_orig_DA = problem_stream.orig_drainage_area
+                max_orig_DA_id = problem_stream.reach_id
+            if problem_stream.orig_drainage_area > problem_stream.fixed_drainage_area:
+                raise Exception("Something weird with the following reach:\n" + str(problem_stream))
+    arcpy.AddMessage("Problem Streams passed check")
+    arcpy.AddMessage("Max problem DA: " + str(max_orig_DA))
+    arcpy.AddMessage("Max problem DA ID: " + str(max_orig_DA_id))
+
+
+def fix_problem_streams(stream_network, problem_streams):
+    """
+    Goes through the stream network and fixes problem streams
+    :param stream_network: The stream network to fix
+    :param problem_streams: A list of problem streams
+    :return:
+    """
+    arcpy.AddMessage("Fixing Streams...")
+    arcpy.AddField_management(stream_network, "Orig_DA", "DOUBLE")
+    req_fields = ["ReachID", "DA_sqkm", "Orig_DA"]
+    with arcpy.da.UpdateCursor(stream_network, req_fields) as cursor:
+        for row in cursor:
+            reach_id = row[0]
+            drain_area = row[1]
+            problem_stream = find_problem_stream(reach_id, problem_streams)
+            if problem_stream:
+                row[1] = problem_stream.fixed_drainage_area
+                row[2] = problem_stream.orig_drainage_area
+            else:
+                row[2] = drain_area
+            cursor.updateRow(row)
+    write_problem_streams(stream_network, problem_streams)
+
+
+def write_problem_streams(stream_network, problem_streams):
+    with open(os.path.join(os.path.dirname(stream_network), "ProblemStreamsList.txt"), 'w') as file:
+        file.write("Something")
+        for problem_stream in problem_streams:
+            file.write("Altered Reach #" + str(problem_stream.reach_id) + '\n')
+
+
+def find_problem_stream(reach_id, problem_streams):
+    """
+    Returns the problem stream that goes with the reach id. If the reach ID is not in the list of problem streams,
+    returns None
+    """
+    for problem_stream in problem_streams:
+        if problem_stream.reach_id == reach_id:
+            return problem_stream
+    return None

RCAT_Stream_Objects.py

@@ -0,0 +1,95 @@
+# -------------------------------------------------------------------------------
+# Name:        StreamObjects
+# Purpose:     This file holds classes that will be used to hold data in RCAT tools
+#
+# Author:      Braden Anderson
+#
+# Created:     03/2018
+# -------------------------------------------------------------------------------
+
+import heapq
+
+
+class DAValueCheckStream:
+    def __init__(self, reach_id,  stream_id, downstream_dist, drainage_area):
+        self.reach_id = reach_id
+        self.stream_id = stream_id
+        self.downstream_dist = downstream_dist
+        self.drainage_area = drainage_area
+
+    def __eq__(self, other):
+        return self.reach_id == other.reach_id
+
+    def __lt__(self, other):
+        """
+        The heap is based on downstream distance, so we define < and > based on downstream distance
+
+        To make the heap a max heap, we reverse what might seem to be intuitive for > and <. To make it a min heap,
+        replace ">" with "<" in the return statement
+        """
+        if not isinstance(other, DAValueCheckStream):
+            raise Exception("Comparing a DAValueCheckStream to another data type is not currently supported")
+        return self.downstream_dist > other.downstream_dist
+
+    def __gt__(self, other):
+        """
+        The heap is based on downstream distance, so we define < and > based on downstream distance
+
+        To make the heap a max heap, we reverse what might seem to be intuitive for > and <. To make it a min heap,
+        replace "<" with ">" in the return statement
+        """
+        if not isinstance(other, DAValueCheckStream):
+            raise Exception("Comparing a DAValueCheckStream to another data type is not currently supported")
+        return self.downstream_dist < other.downstream_dist
+
+    def __str__(self):
+        return str(self.stream_id)
+
+
+class ProblemStream:
+    def __init__(self, reach_id, stream_id, orig_drainage_area, fixed_drainage_area):
+        self.reach_id = reach_id
+        self.stream_id = stream_id
+        self.orig_drainage_area = orig_drainage_area
+        self.fixed_drainage_area = fixed_drainage_area
+
+    def __str__(self):
+        ret_string = 'Reach_ID: ' + str(self.reach_id)
+        ret_string += 'Stream_ID: ' + str(self.stream_id)
+        ret_string += '\nOriginal Drainage Area: ' + str(self.orig_drainage_area)
+        ret_string += '\nFixed Drainage Area: ' + str(self.fixed_drainage_area) + '\n\n'
+
+        return ret_string
+
+
+class StreamHeap:
+    def __init__(self, first_stream):
+        self.streams = [first_stream]
+        self.stream_id = first_stream.stream_id
+
+    def push_stream(self, given_stream):
+        heapq.heappush(self.streams, given_stream)
+
+    def pop(self):
+        return heapq.heappop(self.streams)
+
+    def first_element(self):
+        if len(self.streams) > 0:
+            return self.streams[0]
+        else:
+            return None
+
+    def __eq__(self, other):
+        if not isinstance(other, StreamHeap):
+            raise Exception("A StreamHeap can only be compared to another StreamHeap")
+        return self.stream_id == other.stream_id
+
+    def __str__(self):
+        ret_string = '['
+        for i in range(len(self.streams)):
+            #ret_string += '(' + str(self.streams[i].reach_id) + ', ' + str(self.streams[i].stream_id) + ')'
+            ret_string += str(self.streams[i].downstream_dist)
+            if i + 1 < len(self.streams):
+                ret_string += ', '
+        ret_string += ']'
+        return ret_string

VBET.py

@@ -19,6 +19,7 @@
 import datetime
 import uuid
 import projectxml
+import RCAT_Drainage_Area_Check as DA_Check
 
 
 def main(
@@ -42,7 +43,9 @@ def main(
     scratch,
     ag_distance,
     min_area,
-    min_hole):
+    min_hole,
+    check_drain_area):
+    check_drain_area = parseInputBool(check_drain_area)
 
     arcpy.env.overwriteOutput = True
     arcpy.CheckOutExtension("spatial")
@@ -110,21 +113,21 @@ def main(
     lf = arcpy.ListFields(fcNetwork, "DA_sqkm")
     if len(lf) is 1:
         arcpy.DeleteField_management(fcNetwork, "DA_sqkm")
-    else:
-        pass
 
     arcpy.AddField_management(fcNetwork, "DA_sqkm", "DOUBLE")
-    cursor = arcpy.da.UpdateCursor(fcNetwork, ["MAX", "DA_sqkm"])
-    for row in cursor:
-        row[1] = row[0]
-        cursor.updateRow(row)
-        if row[1] < 0.1:
-            row[1] = 0.1
-        cursor.updateRow(row)
-    del row
-    del cursor
+    with arcpy.da.UpdateCursor(fcNetwork, ["MAX", "DA_sqkm"]) as cursor:
+        for row in cursor:
+            row[1] = row[0]
+            cursor.updateRow(row)
+            if row[1] < 0.1:
+                row[1] = 0.1
+            cursor.updateRow(row)
     arcpy.DeleteField_management(fcNetwork, "MAX")
 
+    if check_drain_area:
+        arcpy.AddMessage("Fixing Drainage Area...")
+        DA_Check.main(fcNetwork)
+
     # create buffers around the different network segments
     arcpy.AddMessage("Creating buffers")
     h = 1
@@ -392,6 +395,13 @@ def main(
     arcpy.CheckInExtension("spatial")
 
 
+def parseInputBool(given_input):
+    if given_input == 'false' or given_input is None:
+        return False
+    else:
+        return True
+
+
 def calc_drain_area(DEM):
     DEMdesc = arcpy.Describe(DEM)
     height = DEMdesc.meanCellHeight