-
Notifications
You must be signed in to change notification settings - Fork 0
/
PNET_Step_4.py
167 lines (121 loc) · 6.45 KB
/
PNET_Step_4.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
import arcpy
import os
from collections import Counter
from PNET_Functions import get_watershed_folders, delete_old, \
finish, delete_temps, remove_fields, attribute_table_to_list, get_field_index
# -------------------------------------------------------------------------------
# Name: PNET Step 4
# Purpose: Edits stream segments and points to create field reaches and points
# Author: Tyler Hatch
#
# Created: 09/23/2019
# Latest Update: 1/22/2020
# -------------------------------------------------------------------------------
# The folder containing all watershed folders
root_folder = arcpy.GetParameterAsText(0)
def main():
# Initialize variables
arcpy.env.overwriteOutput = True
watershed_folders = get_watershed_folders(root_folder)
delete_old(os.path.join(root_folder, "00_ProjectWide", "Intermediates", "Reach_Editing", "Outputs"))
temps_to_delete = []
to_merge_points = []
to_merge_reaches = []
for watershed in watershed_folders:
arcpy.AddMessage("Starting {}...".format(watershed))
# Get file names
input_folder = os.path.join(watershed, "Intermediates", "Reach_Editing", "Inputs")
output_folder = os.path.join(watershed, "Intermediates", "Reach_Editing", "Outputs")
delete_old(output_folder)
stream_seg = os.path.join(input_folder, "Stream_Network_Segments.shp")
points = os.path.join(input_folder, "Points_Merge.shp")
stream_seg_copy = os.path.join(input_folder, "Stream_Network_Segments_To_Edit_Temp.shp")
points_copy = os.path.join(input_folder, "Points_Merge_To_Edit_Temp.shp")
temps_to_delete.append(stream_seg_copy)
temps_to_delete.append(points_copy)
arcpy.Copy_management(stream_seg, stream_seg_copy)
arcpy.Copy_management(points, points_copy)
stream_seg = stream_seg_copy
points = points_copy
# Spatial jon stream network segments by points
fields_to_remove = ["TARGET_FID", "JOIN_FID", "Join_Count"]
spatial_joined = os.path.join(input_folder, "Spatial_Joined_Temp.shp")
temps_to_delete.append(spatial_joined)
remove_fields(fields_to_remove, stream_seg)
arcpy.SpatialJoin_analysis(stream_seg, points, spatial_joined, "JOIN_ONE_TO_MANY")
# Get an attribute table list of the joined shapefile to analyze
stream_seg_list = attribute_table_to_list(spatial_joined)
reach_id_index = get_field_index("TARGET_FID", spatial_joined)
point_id_index = get_field_index("SiteID", spatial_joined)
tor_bor_index = get_field_index("TOR_BOR", spatial_joined)
new_list = split_list_by_id(stream_seg_list, reach_id_index)
to_delete_list = []
keep_points_list = []
# Check which segments we need to delete
for segment in new_list:
if delete_segment(segment, point_id_index, tor_bor_index, reach_id_index):
to_delete_list.append(segment[0][reach_id_index])
else:
keep_points_list.append(segment[0][point_id_index])
segments_layer = "Segments"
arcpy.MakeFeatureLayer_management(stream_seg, segments_layer)
# Delete the segments we need to from initial segments
for to_delete in to_delete_list:
arcpy.SelectLayerByAttribute_management(segments_layer, 'ADD_TO_SELECTION', 'FID = {}'.format(to_delete))
# Save the reaches we want to keep
arcpy.SelectLayerByAttribute_management(segments_layer, 'SWITCH_SELECTION')
reach_save_location = os.path.join(output_folder, "Field_Reaches.shp")
arcpy.CopyFeatures_management(segments_layer, reach_save_location)
to_merge_reaches.append(reach_save_location)
# Save the points we want to keep
points_layer = "Points"
arcpy.MakeFeatureLayer_management(points, points_layer)
for to_keep in keep_points_list:
arcpy.SelectLayerByAttribute_management(points_layer, 'ADD_TO_SELECTION', 'SiteID = {}'.format(to_keep))
point_save_location = os.path.join(output_folder, "Field_Points.shp")
arcpy.CopyFeatures_management(points_layer, point_save_location)
to_merge_points.append(point_save_location)
num_points = int(arcpy.GetCount_management(point_save_location)[0])
num_reaches = int(arcpy.GetCount_management(reach_save_location)[0])
# Check that everything was done correctly
if (num_points/2) != num_reaches:
arcpy.AddMessage("\t This watershed does not have one field reach per two field points!")
arcpy.AddMessage("Saving ProjectWide...")
projectwide_folder = os.path.join(root_folder, "00_ProjectWide", "Intermediates", "Reach_Editing", "Outputs")
arcpy.Merge_management(to_merge_points, os.path.join(projectwide_folder, "Field_Points.shp"))
arcpy.Merge_management(to_merge_reaches, os.path.join(projectwide_folder, "Field_Reaches.shp"))
delete_temps(temps_to_delete)
finish()
def split_list_by_id(segment_list, id_index):
to_return = []
current_seg = []
previous_id = segment_list[0][id_index]
for segment in segment_list:
# Adding a segment to the same reach
if segment[id_index] == previous_id:
current_seg.append(segment)
# Adding a new reach to the return list, then creating a new reach
else:
to_return.append(current_seg)
current_seg = [segment]
previous_id = segment[id_index]
# This makes sure that the final reach gets added, because it may be missed in the for loop
to_return.append(current_seg)
# Returns a list of reaches, which each contain a list of segments, which each contain a list of field values
return to_return
def delete_segment(segment, point_id, tor_bor, reach_id):
# If this segment is not touching exactly two points, delete it
if len(segment) != 2:
return True
# These hold data for the first and second point that the segment is touching
point_a = segment[0]
point_b = segment[1]
# If the two points the segment is touching are from different sites, delete it
if point_a[point_id] != point_b[point_id]:
return True
# Also delete if they are both TOR or both BOR
if point_a[tor_bor] == point_b[tor_bor]:
return True
return False
if __name__ == "__main__":
main()