|
|
楼主 |
发表于 2018-12-26 20:14
|
显示全部楼层
工具代码如下：
import sys, arcpy, traceback, math, string, os, fileinput, numpy, operator, scipy
from scipy.spatial import cKDTree as KDTree
from scipy import ndimage
from operator import itemgetter
from arcpy.sa import *
from arcpy import env
import arcpy.cartography as CA
arcpy.env.overwriteOutput = True
arcpy.Delete_management("in_memory")
glacier_outlines = arcpy.GetParameterAsText(0)
dem = arcpy.GetParameterAsText(1)
tributary_threshold1 = int(arcpy.GetParameterAsText(2))
tributary_threshold2 = int(arcpy.GetParameterAsText(3))
flow_line = arcpy.GetParameterAsText(4)
cellsize = arcpy.GetRasterProperties_management(dem,"CELLSIZEX")
cellsize_int = int(cellsize.getOutput(0))
arcpy.env.outputCoordinateSystem = dem
path_split_flow_line = os.path.split(flow_line)
flow_line_path = path_split_flow_line[0]
flow_line_name = path_split_flow_line[1]
outline_count = arcpy.GetCount_management(glacier_outlines) #COUNT HOW MANY OUTLINES IN TOTAL (FOR PROGRESS REPORT)
outline_count_result = int(outline_count.getOutput(0)) #GET INT RESULT OF OUTLINE COUNT FOR PROGRESS REPORT)
progress_counter = 1
def create_axis(outline):
outline_feature_count = arcpy.GetCount_management(outline) #check outline is a single feature - count"
outline_feature_count_result = int(outline_feature_count.getOutput(0)) #check outline is a single feature - get result"
if outline_feature_count_result != 1: #check outline is a single feature"
outline = arcpy.Dissolve_management(outline, "in_memory\\dissolved_outline")
mbg = arcpy.MinimumBoundingGeometry_management(outline, "in_memory\\mbg", "CONVEX_HULL", "NONE", "","MBG_FIELDS") #create minimum bounding geometry, convex hull method"
axis = arcpy.XYToLine_management(mbg, "in_memory\\axis_out","MBG_APodX1", "MBG_APodY1", "MBG_APodX2", "MBG_APodY2") # Create long axis from fields in mbg
return axis
def extend_line(original_line, extend_percentage):
arcpy.env.extent = "MAXOF"
with arcpy.da.SearchCursor(original_line, ["SHAPE@", "SHAPE@LENGTH"]) as cursor: #set up search cursor to access geometry of line
for row in cursor:
extend_length_axis = row[1]*extend_percentage
firstpoint = row[0].firstPoint #access firstpoint co-ords
lastpoint = row[0].lastPoint #access lastpoint co-ords
if firstpoint.X - lastpoint.X == 0: #test to see if line is vertical (infinite slope so calculating slope would result in error)
new_firstpoint_X = firstpoint.X #Line is vertical so new x co-ordinates will be the same as original
new_lastpoint_X = lastpoint.X #Line is vertical so new x co-ordinates will be the same as original
if firstpoint.Y > lastpoint.Y: #Test to see if firstpoint above lastpoint
new_firstpoint_Y = firstpoint.Y + extend_length_axis #Firstpoint is above, so add the extend length
new_lastpoint_Y = lastpoint.Y - extend_length_axis #Lastpoint is below, so minus the extend length
else: #test if firstpoint is below
new_firstpoint_Y = firstpoint.Y - extend_length_axis #firstpoint is below, so minus the extend length
new_lastpoint_Y = lastpoint.Y + extend_length_axis #lastpoint is above so add extend length
else: #Line is not vertical, so slope can be calculated
axis_slope = ((firstpoint.Y - lastpoint.Y)/(firstpoint.X - lastpoint.X)) #Calculate slope (change in y/change in x)
if axis_slope <0: #Test if slope is negative
axis_slope = axis_slope*-1 #Invert negative slopes to make calculations simpler
axis_slope_radians = math.atan(axis_slope) #Convert slope to radians
changex = float(math.cos(axis_slope_radians))*float(extend_length_axis) #Calculate amount to change X co-ord by (cos of slope * extend length)
changey = float(math.sin(axis_slope_radians))*float(extend_length_axis) #Calculate amount to change Y co-ord by (sin of slope * extend length)
if firstpoint.X > lastpoint.X: #Test if firstpoint X co-ord is greater than lastpoint (ie to the right)
new_firstpoint_X = firstpoint.X + changex #Firstpoint X co-ord is greater (to the right) of lastpoint, so add changeX
new_lastpoint_X = lastpoint.X - changex #Firstpoint X co-ord is smaller (to the left) of lastpoint, so minus changeX
else: #Test if firstpoint X co-ord is smaller (left) of lastpoint
new_firstpoint_X = firstpoint.X - changex #Firstpoint X co-ord is smaller (left) of lastpoint, so minus changex
new_lastpoint_X = lastpoint.X + changex #Lastpoint X co-ord is greater (right) of firstpoint, so add changex
if firstpoint.Y >= lastpoint.Y: #Test if firstpoint X co-ord is greater (above) lastpoint. Equal implies a horizontal line
new_firstpoint_Y = firstpoint.Y + changey #Firstpoint is above, so add change Y
new_lastpoint_Y = lastpoint.Y - changey #Lastpoint is below so minus change Y
else: #Test if lastpoint is above firstpoint
new_firstpoint_Y = firstpoint.Y - changey #firstpoint is below, so minus change Y
new_lastpoint_Y = lastpoint.Y + changey #lastpoint is above so add change Y
pointfc = arcpy.CreateFeatureclass_management("in_memory","extend_line_pt","POINT","","","",original_line) #create new feature class (in memory) for new xy points
point_cursor = arcpy.da.InsertCursor(pointfc, ["SHAPE@XY"]) #set up insert cursor to populate feature class
new_firstpoint = (new_firstpoint_X, new_firstpoint_Y) #gather firstpoint co-ords
new_lastpoint =(new_lastpoint_X, new_lastpoint_Y) #gather lastpoint co-ords
point_cursor.insertRow([new_firstpoint]) #populate firstpoint
point_cursor.insertRow([new_lastpoint]) #populate lastpoint
extended_line = arcpy.PointsToLine_management(pointfc, "in_memory\\extended_line_line") #draw line (new extended axis) between first and last points
return extended_line
def points_on_line(line, resolution):
new_points = arcpy.CreateFeatureclass_management("in_memory", "points","POINT")
new_points_cursor = arcpy.da.InsertCursor(new_points, ('SHAPE@'))
start = 0
with arcpy.da.SearchCursor(line, ["SHAPE@LENGTH", "SHAPE@"]) as cursor:
for row in cursor:
length = row[0]
while start < length:
new_point = row[1].positionAlongLine(start)
new_points_cursor.insertRow((new_point,))
start = start + resolution
del new_points_cursor
return new_points
def create_perpendiculars(in_lines, distance):
perpendiculars = arcpy.CreateFeatureclass_management("in_memory", "perps","POLYLINE","","","",dem)
new_line_cursor = arcpy.da.InsertCursor(perpendiculars, ('SHAPE@'))
with arcpy.da.SearchCursor(in_lines, ["SHAPE@"]) as cursor:
for row in cursor:
startx = row[0].firstPoint.X
starty = row[0].firstPoint.Y
endx = row[0].lastPoint.X
endy = row[0].lastPoint.Y
midx = row[0].centroid.X
midy = row[0].centroid.Y
if starty==endy or startx==endx:
if starty == endy:
y1 = starty + distance
y2 = starty - distance
x1 = startx
x2 = startx
if startx == endx:
y1 = starty
y2 = starty
x1 = startx + distance
x2 = startx - distance
else:
m = ((starty - endy)/(startx - endx)) #get the slope of the line
negativereciprocal = -1*((startx - endx)/(starty - endy)) #get the negative reciprocal
if m > 0:
if m >= 1:
y1 = negativereciprocal*(distance)+ starty
y2 = negativereciprocal*(-distance) + starty
x1 = startx + distance
x2 = startx - distance
if m < 1:
y1 = starty + distance
y2 = starty - distance
x1 = (distance/negativereciprocal) + startx
x2 = (-distance/negativereciprocal)+ startx
if m < 0:
if m >= -1:
y1 = starty + distance
y2 = starty - distance
x1 = (distance/negativereciprocal) + startx
x2 = (-distance/negativereciprocal)+ startx
if m < -1:
y1 = negativereciprocal*(distance)+ starty
y2 = negativereciprocal*(-distance) + starty
x1 = startx + distance
x2 = startx - distance
array = arcpy.Array([arcpy.Point(x1,y1),arcpy.Point(x2, y2)])
polyline = arcpy.Polyline(array)
new_line_cursor.insertRow([polyline])
return perpendiculars
def central_points_on_perps(perps, glacier_outline, dem):
clipped_axis_perps = arcpy.Clip_analysis(perps, glacier_outline, "in_memory\\clipped_axis_perps")
split_axis_perp_lines = arcpy.SplitLine_management(clipped_axis_perps, "in_memory\\split_axis_perps")
central_points_no_alt = arcpy.FeatureToPoint_management (split_axis_perp_lines, "in_memory\\central_no_alt")
central_points_with_alt = ExtractValuesToPoints(central_points_no_alt, dem, "in_memory\\central_points_with_alt")
return central_points_with_alt
def flowline (central_points_with_alt, flow_line_output, flowline_glacier_outline, smooth_tolerance):
inaccessible_lowpoint = 0
glacier_outline_line = arcpy.FeatureToLine_management(flowline_glacier_outline, "in_memory\\glacier_outline_line")
dissolve_outline_line = arcpy.Dissolve_management(glacier_outline_line, "in_memory\\dissolve_glacier_outline_line")
outline_line_geom = arcpy.CopyFeatures_management(dissolve_outline_line, arcpy.Geometry())
fields = ["RASTERVALU", "SHAPE@XY"] #list of fields used in central points: row[0] = elevation, row[1] = co-ords
central_point_dict = {} #create empty dictionary for all central points
flowline_points_dict = {}
max_min_list = []
with arcpy.da.SearchCursor(central_points_with_alt, "RASTERVALU") as cursor:
for row in cursor:
max_min_list.append(row[0])
max_central_point_remove = max(max_min_list)
min_central_point_remove = min(max_min_list)
with arcpy.da.UpdateCursor(central_points_with_alt, "RASTERVALU") as cursor: ##DELETE MAX AND MIN POUNTS TO STOP CODE FAILING CLOSE TO MARGIN
for row in cursor:
if row[0] == max_central_point_remove:
cursor.deleteRow()
if row[0] == min_central_point_remove:
cursor.deleteRow()
with arcpy.da.SearchCursor(central_points_with_alt, fields) as cursor:
for row in cursor:
if row[0] != -9999: #Excludes any 'no data' points
index = row[1] # use elevation as index
value = row[0] # use co-ords as values
central_point_dict[row[1]] = row[0] #populate dictionary
maximum_central_point_alt = max(central_point_dict.itervalues()) #get max value
min_central_point_alt = min(central_point_dict.itervalues()) #get min value
max_co_ords_list = [] #set up empty list for high point(s)
for index, value in central_point_dict.iteritems():
if value == maximum_central_point_alt:
max_co_ords_list.append(index)
min_co_ords_list = [] #set up empty list for high point(s)
for index, value in central_point_dict.iteritems():
if value == min_central_point_alt:
min_co_ords_list.append(index)
end_co_ords = min_co_ords_list[0]
if len(max_co_ords_list) > 1: #test if more than 1 high point
arcpy.AddMessage("WARNING multiple high points, using first point")
start_ele = maximum_central_point_alt #set initial starting elevation to maximum
start_co_ords = max_co_ords_list[0] #set initial starting co-ords to that of maximum
flowline_points_dict = {} # Dictionary to hold initial flowline points in
line_segment_list = [] # list to hold all line segments in
lowpoint_stop = 0 #boolean to contineue or stop
while (lowpoint_stop == 0):
with arcpy.da.SearchCursor(central_points_with_alt, fields) as cursor:
nearest_neighbour_list = [] #create new empty list for each loop
for row in cursor:
if row[0] < start_ele and row[0] != -9999: #find all lower points (excluding nodata)
nearest_neighbour_list.append(row[1]) #make list of x co ord and y co ord
if start_co_ords not in nearest_neighbour_list:
nearest_neighbour_list.append(start_co_ords) #add starting point to list
len_nearest_neighbour_list = len(nearest_neighbour_list) #get length of list
if len_nearest_neighbour_list == 0:
lowpoint_stop = 1 #stop as no more lower nearest neighbours
else:
index_val = nearest_neighbour_list.index(start_co_ords)
ArrayOfPoints = numpy.array(nearest_neighbour_list) #convert list to numpy array
KDTreeOfPoints = KDTree(ArrayOfPoints)
closest_neighbour = 2 #closest neighbour = 2 as 1 = original)
return_neighbour = 1 #closets neighbour
distances, indicies = KDTreeOfPoints.query(ArrayOfPoints[index_val],closest_neighbour) #closest neighbour = 2 as 1 = original)
NearestNeighbor = (nearest_neighbour_list[indicies[return_neighbour]])
pointlist = [] #create new list to hold start and end point for segment
pointlist.append([start_co_ords,NearestNeighbor]) #append start and end point
features = []
for feature in pointlist:
features.append(arcpy.Polyline(arcpy.Array([arcpy.Point(*coords) for coords in feature]),flowline_glacier_outline))
for segment in features:
intersect = segment.crosses(outline_line_geom[0])
if intersect == 0:
line_segment_list.append(segment)
else:
while intersect == 1:
if NearestNeighbor != end_co_ords:
closest_neighbour = closest_neighbour + 1
return_neighbour = return_neighbour + 1
distances, indicies = KDTreeOfPoints.query(ArrayOfPoints[index_val],closest_neighbour) #closest neighbour = 2 as 1 = original)
NearestNeighbor = (nearest_neighbour_list[indicies[return_neighbour]])
pointlist_intersect = [] #create new list to hold start and end point for segment
pointlist_intersect.append([start_co_ords,NearestNeighbor]) #append start and end point
features_intersect = []
for feature in pointlist_intersect:
features_intersect.append(arcpy.Polyline(arcpy.Array([arcpy.Point(*coords) for coords in feature]),flowline_glacier_outline))
for segment_intersect in features_intersect:
intersect = segment_intersect.crosses(outline_line_geom[0])
if intersect == 0:
line_segment_list.append(segment_intersect)
else:
intersect = 0
start_co_ords = NearestNeighbor #reset start co-ords to that of new nearest neighbour
start_ele = central_point_dict[NearestNeighbor] #reset start elevation to that of new nearest neighbour
closest_neighbour = 2 #closest neighbour = 2 as 1 = original)
return_neighbour = 1 #closest neighbour
initial_flowline = arcpy.CopyFeatures_management(line_segment_list, "in_memory\\initial_flow_line")
dissolved_flowline = arcpy.Dissolve_management(initial_flowline, "in_memory\\dissolved_flowline")
smooth_flowline = CA.SmoothLine(dissolved_flowline, "in_memory\\flow_line", "PAEK", smooth_tolerance)
smooth_flowline_geometry = arcpy.CopyFeatures_management(smooth_flowline, arcpy.Geometry())
for geometry in smooth_flowline_geometry:
intersect = geometry.crosses(outline_line_geom[0])
if intersect == 1:
arcpy.AddMessage("WARNING: SMOOTHING CAUSED CENTRELINE TO CROSS OUTLINE")
return smooth_flowline
def new_branch(input_flowline, branch_outline, input_outline):
clipped_dem = arcpy.Clip_management(dem,"","in_memory\\clipped_dem",input_outline,"","ClippingGeometry")
clipped_dem_raster = Raster(clipped_dem)
up_list = []
up_list = []
arcpy.env.extent = clipped_dem_raster
global new_branch_count
new_branch_count = 0
dem_numpy = arcpy.RasterToNumPyArray(clipped_dem_raster,"","","",0)
d8_structure = ndimage.generate_binary_structure(2, 2)
upstream_area_original = 1
glacier_numpy = dem_numpy.copy()
glacier_numpy[glacier_numpy != 0] = 1
with arcpy.da.SearchCursor(input_flowline, ["SHAPE@LENGTH"]) as cursor:
for row in cursor:
flowline_length = row[0]
flowline_point_res = flowline_length/100.0
flowline_points = points_on_line(input_flowline, flowline_point_res)
branch_points_alt = ExtractValuesToPoints(flowline_points, dem, "in_memory\\branch_points_altitude")
branch_points_alt_layer = arcpy.MakeFeatureLayer_management(branch_points_alt, "in_memory\\branch_points_layer")
with arcpy.da.SearchCursor(branch_points_alt_layer, [arcpy.Describe(branch_points_alt_layer).OIDFieldName, "RASTERVALU"]) as cursor:
for row in cursor:
if row[0] > 1:
if new_branch_count == 0:
query = arcpy.Describe(branch_points_alt_layer).OIDFieldName+"= "+str(row[0])
arcpy.SelectLayerByAttribute_management(branch_points_alt_layer, "NEW_SELECTION", query)
raster_flowline_general = Raster(arcpy.PointToRaster_conversion(branch_points_alt_layer, "RASTERVALU", "in_memory\\ras1", "","", cellsize_int))
flowline_numpy = arcpy.RasterToNumPyArray(raster_flowline_general,"","","",0)
original_cell = numpy.amax(flowline_numpy)
itemindex = numpy.where(flowline_numpy == original_cell)
row = itemindex[0]
col = itemindex[1]
dem_copy = dem_numpy.copy()
dem_copy[dem_copy < original_cell] = 0
dem_copy[dem_copy >= original_cell] = 1
labeled_array, num_features = scipy.ndimage.measurements.label(dem_copy, structure=d8_structure)
zone = labeled_array[row,col]
zone_cells = numpy.argwhere(labeled_array == zone)
total_cells = int(len(zone_cells))
upstream_area = total_cells*cellsize_int*cellsize_int
upstream_area_change = upstream_area - upstream_area_original
percentage_change = (float(upstream_area)/float(upstream_area_original))*100
up_list.append(upstream_area)
upstream_area_original = upstream_area
if upstream_area < area_threshold2 or percentage_change < (100.0 + tributary_threshold1) or percentage_change >500:
newarray = labeled_array.copy()
newarray[newarray != zone] = 0
newarray[newarray == zone] = 1
glacier_numpy[newarray == 1] = 0
else:
search_branch = 1
left_co_ord = arcpy.GetRasterProperties_management(clipped_dem,"LEFT")
bottom_co_ord = arcpy.GetRasterProperties_management(clipped_dem,"BOTTOM")
left_co_ord_float = float(left_co_ord.getOutput(0))
bottom_co_ord_float = float(bottom_co_ord.getOutput(0))
bottom_left_point = arcpy.Point(left_co_ord_float,bottom_co_ord_float)
new_branch_raster = Int(arcpy.NumPyArrayToRaster(glacier_numpy, bottom_left_point, cellsize_int, cellsize_int, 0))
new_branch_outline_unclipped = arcpy.RasterToPolygon_conversion(new_branch_raster, "in_memory\\new_branch_outline_unclipped", "SIMPLIFY")
new_branch_outline = arcpy.Clip_analysis(new_branch_outline_unclipped, glacier_outline, branch_outline)
new_branch_count = 1
return new_branch_outline
flow_line_final = arcpy.CreateFeatureclass_management(flow_line_path, flow_line_name, "POLYLINE")
with arcpy.da.SearchCursor(glacier_outlines, [arcpy.Describe(glacier_outlines).OIDFieldName]) as cursor:
for row in cursor:
arcpy.env.extent = glacier_outlines
glacier_outline_layer = arcpy.MakeFeatureLayer_management(glacier_outlines, "in_memory\\glacier_outlines_layer")
flow_line_final_individual = arcpy.CreateFeatureclass_management("in_memory","flow_line_final_individual", "POLYLINE")
arcpy.AddMessage("Generating centreline(s) for glacier "+str(progress_counter)+" of "+str(outline_count_result))
query = arcpy.Describe(glacier_outline_layer).OIDFieldName+"="+str(row[0])
arcpy.SelectLayerByAttribute_management(glacier_outline_layer, "NEW_SELECTION", query)
glacier_outline = arcpy.CopyFeatures_management(glacier_outline_layer, "in_memory\\layer_test_hope")
with arcpy.da.SearchCursor(glacier_outline, ["SHAPE@AREA"]) as cursor:
for row in cursor:
outline_area = row[0]
area_threshold2 = outline_area * (tributary_threshold2 / 100.0)
axis = create_axis(glacier_outline)
extended_axis = extend_line(axis, 0.25)
points_on_axis = points_on_line(extended_axis, cellsize_int*3)
split_axis = arcpy.SplitLineAtPoint_management (extended_axis, points_on_axis, "in_memory\\splitaxis", 1)
perpendiculars = create_perpendiculars(split_axis, 1000000)
central_points_with_alt = central_points_on_perps(perpendiculars, glacier_outline, dem)
initial_flowline = flowline(central_points_with_alt, "in_memory\\initial_flow_line", glacier_outline, 1000)
new_branch_poly = new_branch(initial_flowline, "in_memory\\new_outline", glacier_outline)
if new_branch_count == 0:
arcpy.Append_management(initial_flowline, flow_line_final)
while new_branch_count == 1:
arcpy.AddMessage("Multiple Branches Detected")
secondary_axis = create_axis(new_branch_poly)
secondary_extended_axis = extend_line(secondary_axis, 0.25)
secondary_points_on_axis = points_on_line(secondary_extended_axis, cellsize_int*3)
secondary_split_axis = arcpy.SplitLineAtPoint_management (extended_axis, secondary_points_on_axis, "in_memory\\splitaxis", 1)
secondary_perpendiculars = create_perpendiculars(secondary_split_axis, 1000000)
secondary_central_points_with_alt = central_points_on_perps(secondary_perpendiculars, new_branch_poly, dem)
secondary_flowline = flowline(secondary_central_points_with_alt, "in_memory\\secondary_flow_line", new_branch_poly, 1000)
secondary_flowline_copied = arcpy.CopyFeatures_management(secondary_flowline, "in_memory\\secondary_flow_line_copied")
corrected_poly = new_branch(secondary_flowline, "in_memory\\corrected_poly1", glacier_outline)
corrected_axis = create_axis(corrected_poly)
corrected_extended_axis = extend_line(corrected_axis, 0.25)
corrected_points_on_axis = points_on_line(corrected_extended_axis, cellsize_int*3)
corrected_split_axis = arcpy.SplitLineAtPoint_management (corrected_extended_axis, corrected_points_on_axis, "in_memory\\corrected_splitaxis", 1)
corrected_perpendiculars = create_perpendiculars(corrected_split_axis, 1000000)
corrected_central_points_with_alt = central_points_on_perps(corrected_perpendiculars, corrected_poly, dem)
corrected_flowline = flowline(corrected_central_points_with_alt, "in_memory\\corrected_additional_flowline", corrected_poly, 1000)
corrected_flowline_copied = arcpy.CopyFeatures_management(corrected_flowline, "in_memory\\corrected_flowline_copied")
flowline_intersect_points = arcpy.Intersect_analysis([corrected_flowline_copied, secondary_flowline_copied],"in_memory\\flow_intersect_points","","","POINT")
secondary_fl_split = arcpy.SplitLineAtPoint_management(secondary_flowline_copied, flowline_intersect_points, "in_memory\\secondary_fl_split", 5)
secondary_fl_split_layer = arcpy.MakeFeatureLayer_management(secondary_fl_split, "in_memory\\secondary_fl_split_layer")
max_midpoint = 0.0
with arcpy.da.UpdateCursor(secondary_fl_split_layer, arcpy.Describe(secondary_fl_split_layer).OIDFieldName) as cursor:
for row in cursor:
query = arcpy.Describe(secondary_fl_split_layer).OIDFieldName+"="+str(row[0])
arcpy.SelectLayerByAttribute_management(secondary_fl_split_layer, "NEW_SELECTION", query)
midpoint_seecondary_fl_split = arcpy.FeatureVerticesToPoints_management(secondary_fl_split_layer, "in_memory\\midpoint_seecondary_fl_split", "MID")
midpoint_seecondary_fl_split_ele = ExtractValuesToPoints(midpoint_seecondary_fl_split, dem, "in_memory\\midpoint_seecondary_fl_split_elevation")
with arcpy.da.UpdateCursor(midpoint_seecondary_fl_split_ele, "RASTERVALU") as point_cursor:
for point_row in point_cursor:
elevation = point_row[0]
if elevation > max_midpoint:
max_midpoint = elevation
else:
cursor.deleteRow()
arcpy.SelectLayerByAttribute_management(secondary_fl_split_layer, "CLEAR_SELECTION")
new_branch_poly = new_branch(secondary_flowline, new_branch_poly, new_branch_poly)
arcpy.Append_management([secondary_fl_split_layer,corrected_flowline_copied], flow_line_final_individual)
centreline_count_glacier = arcpy.GetCount_management(flow_line_final_individual) #COUNT HOW MANY centrelines IN glacier
centreline_count_result = int(centreline_count_glacier.getOutput(0)) #GET INT RESULT
split_points = arcpy.FeatureVerticesToPoints_management(flow_line_final_individual, "in_memory\\flowline_split_points", "BOTH_ENDS")
split_flolwine_glacier_ori = arcpy.SplitLineAtPoint_management (flow_line_final_individual, split_points, "in_memory\\split_flolwine_glacier_ori", 1)
split_flolwine_glacier = arcpy.MakeFeatureLayer_management(split_flolwine_glacier_ori, "in_memory\\split_flolwine_glacier")
start_ele_list = []
start_ele_dict = {}
with arcpy.da.SearchCursor(split_flolwine_glacier, [arcpy.Describe(split_flolwine_glacier).OIDFieldName]) as cursor:
for row in cursor:
flowline_id = row[0]
flowline_query = arcpy.Describe(split_flolwine_glacier).OIDFieldName+" = "+str(flowline_id)
arcpy.SelectLayerByAttribute_management(split_flolwine_glacier,"NEW_SELECTION",flowline_query)
line_startpoint = arcpy.FeatureVerticesToPoints_management(split_flolwine_glacier, "in_memory\\line_in_startpoint", "START")
line_endpoint = arcpy.FeatureVerticesToPoints_management(split_flolwine_glacier, "in_memory\\line_in_endpoint", "END")
startpoint_ele_point = ExtractValuesToPoints(line_startpoint, dem,"in_memory\\startpoint_point_ele") #Get elevation of each point
endpoint_ele_point = ExtractValuesToPoints(line_endpoint, dem,"in_memory\\endpoint_point_ele") #Get elevation of each point
with arcpy.da.SearchCursor(startpoint_ele_point, "RASTERVALU") as cursor2:
for row2 in cursor2:
startpoint_ele = row2[0]
with arcpy.da.SearchCursor(endpoint_ele_point, "RASTERVALU") as cursor3:
for row3 in cursor3:
endpoint_ele = row3[0]
if startpoint_ele > endpoint_ele:
arcpy.FlipLine_edit(split_flolwine_glacier)
startpoint_checked = line_endpoint
endpoint_checked = line_startpoint
startpoint_ele_checked = endpoint_ele
endpoint_ele_checked = startpoint_ele
else:
startpoint_checked = line_startpoint
endpoint_checked = line_endpoint
startpoint_ele_checked = startpoint_ele
endpoint_ele_checked = endpoint_ele
start_ele_list.append(startpoint_ele_checked)
start_ele_dict[row[0]] = startpoint_ele_checked
min_startpoint_ele = min(start_ele_list)
arcpy.SelectLayerByAttribute_management(split_flolwine_glacier,"CLEAR_SELECTION")
dangle_points = arcpy.FeatureVerticesToPoints_management(split_flolwine_glacier, "in_memory\\dangle_points", "DANGLE")
start_points = arcpy.FeatureVerticesToPoints_management(split_flolwine_glacier, "in_memory\\start_points", "START")
dangle_points_layer = arcpy.MakeFeatureLayer_management(dangle_points, "in_memory\\dangle_points_layer")
start_points_layer = arcpy.MakeFeatureLayer_management(start_points, "in_memory\\start_points_layer")
arcpy.SelectLayerByLocation_management(dangle_points_layer, "ARE_IDENTICAL_TO", start_points_layer)
arcpy.SelectLayerByLocation_management(dangle_points_layer, "","","","SWITCH_SELECTION")
dangle_count_result = arcpy.GetCount_management(dangle_points_layer)
dangle_count = int(dangle_count_result.getOutput(0))
arcpy.SelectLayerByLocation_management(split_flolwine_glacier, "BOUNDARY_TOUCHES", dangle_points_layer)
with arcpy.da.SearchCursor(split_flolwine_glacier, [arcpy.Describe(split_flolwine_glacier).OIDFieldName]) as cursor:
for row in cursor:
start_ele = start_ele_dict[row[0]]
dangle_query = arcpy.Describe(split_flolwine_glacier).OIDFieldName+" = "+str(row[0])
arcpy.SelectLayerByAttribute_management(split_flolwine_glacier,"NEW_SELECTION",dangle_query)
start_segment = arcpy.CopyFeatures_management(split_flolwine_glacier, "in_memory\\single_dangle")
single_flowline = arcpy.CopyFeatures_management(split_flolwine_glacier, "in_memory\\flowline")
while start_ele > min_startpoint_ele:
arcpy.SelectLayerByLocation_management(split_flolwine_glacier, "INTERSECT", start_segment)
with arcpy.da.SearchCursor(split_flolwine_glacier, [arcpy.Describe(split_flolwine_glacier).OIDFieldName]) as cursor:
for row in cursor:
seg_ele = start_ele_dict[row[0]]
if seg_ele < start_ele:
seg_select_query = arcpy.Describe(split_flolwine_glacier).OIDFieldName+" ="+str(row[0])
arcpy.SelectLayerByAttribute_management(split_flolwine_glacier,"NEW_SELECTION",seg_select_query)
arcpy.Append_management(split_flolwine_glacier, single_flowline)
start_segment = arcpy.CopyFeatures_management(split_flolwine_glacier, "in_memory\\new_start")
start_ele = seg_ele
single_flowline_dissolve = arcpy.Dissolve_management(single_flowline, "in_memory\\single_flowline_dissolve")
arcpy.Append_management(single_flowline_dissolve, flow_line_final)
progress_counter = progress_counter + 1
arcpy.Delete_management("in_memory")
|
|
雷达卡
发表于 2018-12-26 20:13
提升卡
沉默卡
喧嚣卡
变色卡
抢沙发
显身卡
发表于 2018-12-27 09:04
