-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathRCAT_Drainage_Area_Check.py
170 lines (144 loc) · 6.82 KB
/
RCAT_Drainage_Area_Check.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
# -------------------------------------------------------------------------------
# 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:
"""
arcpy.AddMessage("Fixing drainage area...")
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 streams that need drainage area updated...")
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("Updating drainage area values...")
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