PrivacyScanData/truncate/new2/Polygon-Dilation-master/corners_narrows_workflow.py

# -*- coding: utf-8 -*-
"""
This is a workflow to truncate deep and narrow corners and dilate narrow segments 
of soil polygons. This workflow is intened to assist in the QA/QC process of SSURGO

Five genveral steps:
    1) Smooth and Generalize features
    2) Identify deep corners and short narrow ends
    3) Truncate features which only share boundary with two MU's
    4) Identify narrow segements of polygons
    5) Dilate narrow segments and update

@author: alexander.stum

Last revised 5/30/18

PrintMsg and errorMsg were written by Steve Peaslee. At the momment these functions
are not utilized. When set-up to be used as an ArcToolbox Script, they will become
useful


"""

import arcpy,traceback,sys

#======= User Defined Variables  ==========
#The geodatabase with MUPOLYGON and SAPOLYGON
gdb = r'G:\GIS_specialist\Data\Geospatial\soils\FY18\Lubbock\Borden\workspace\Lubbock18_aks051618.gdb'
arcpy.env.workspace = gdb
#The main soil polygon layer
mupolygon = "MUPOLYGON"
#The survey area boundary, NOTE that no features touching boundary will be modified
sapolygon = "SAPOLYGON"
#selection query, specify if the entire survey area is not being evaluated
#otherwise set to None
query = "Update_Name IS NOT NULL"
#query = None

#======= Set Parameters  ==========
#Parameter for line smooting of polylines. Subjective, I like 25
smooth_tol = "20 Meters"
#Parameter for the generalization of polylines. Subjective, I like 2.
gen_tol = "2 Meters"
#Half the miniumum width standard, currently national standard is 38 m
min_width = 19
#All deep/narrow corners with Shape_Length less than this parameter will not 
#be treated. Leave it to the human
min_len = "300"
#Look ahead parameter, used by centerline function
f = 5


#======= Functions  ==========
def PrintMsg(msg, severity=0):
    # Adds tool message to the geoprocessor and prints to console window

    #Split the message on \n first, so that if it's multiple lines, a GPMessage will be added for each line
    try:
        for string in msg.split('\n'):
            #Add a geoprocessing message (in case this is run as a tool)
            if severity == 0:
                arcpy.AddMessage(string)

            elif severity == 1:
                arcpy.AddWarning(string)

            elif severity == 2:
                arcpy.AddMessage("    ")
                arcpy.AddError(string)

    except:
        pass


## ===================================================================================
def errorMsg():
    try:
        tb = sys.exc_info()[2]
        tbinfo = traceback.format_tb(tb)[0]
        theMsg = tbinfo + "\n" + str(sys.exc_type)+ ": " + str(sys.exc_value)
        PrintMsg(theMsg, 2)

    except:
        PrintMsg("Unhandled error in errorMsg method", 2)
        pass
    
    
## ===================================================================================   
def centerline(narrows, centerline, extend_to, f=5):
    """Calculates centerline between two more or less parrallel polylines (multipart)
    which represent the outline of a polygon.
    
    It finds the centerline by starting at the endpoint of a feature. It first
    looks at the first five vertices of the opposing part. If the nearest vertex
    is the last of these five, it will look at the next five. It does this for
    each vertex of the first part. At the end it verifies that the enddpoints of
    the second part were involved, if not it similarly finds the nearest vertices
    in part one to the uninvolved endbpoints of part 2.
    As the two part polyline outline is produced by buffering, the algorithm extends
    the centerline to the nearest negative buffer feature to extend the centerline
    beyond the inward bowing minicus.

    - **Parameters**::
    
        Name:      Type:           Description:
        narrows     String          Name of the polyline outline feature
        centerline  String          Name of the centerline feature created by function
        extend_to   String          Name of the polygon feature to which centerline will be stretched to
        f           Integer         Look ahead parameter, how many verticies to consider at a time
    
    - **Returns**::
        
        True        if successful
        
    .. note:: Requires ArcPy
`
    """  
    
    s_r = arcpy.Describe(narrows).spatialReference                  
    arcpy.CreateFeatureclass_management(arcpy.env.workspace, centerline, 'POLYLINE',\
                                        spatial_reference= s_r)
    
    def midpoint(pt1,pt2):
        x = pt1.X-(pt1.X-pt2.X)/2.0
        y = pt1.Y-(pt1.Y-pt2.Y)/2.0
        return arcpy.Point(x,y)
    
    sCur = arcpy.da.SearchCursor(narrows,'SHAPE@')
    iCur = arcpy.da.InsertCursor(centerline,'SHAPE@')
    for poly, in sCur:
        if poly and poly.partCount==2:
            line  = arcpy.Array()
            pi = 0   #current position in part 2
            p1 = poly.getPart(0)
            p2 = [p for p in poly.getPart(1)] #put in list to be able to slice
            p1c = p1.count
            dlast0 = ((p1[p1c-1].X-p2[0].X)**2+(p1[p1c-1].Y-p2[0].Y)**2)**.5
            d00= ((p1[0].X-p2[0].X)**2+(p1[0].Y-p2[0].Y)**2)**.5
            if dlast0 < d00:
                index=  range(p1c-1,-1,-1)
                line.add(midpoint(p2[0],p1[p1.count-1]))
            else:
                index=  range(p1c)
                line.add(midpoint(p2[0],p1[0]))           
    
            for i in index:
                pf = f-1 #relative forward position
                while pf == f-1:        #if the relative forward position is the same as the last relative position, continue looking
                    d = [((p1[i].X-p.X)**2+(p1[i].Y-p.Y)**2)**.5 \
                         for p in p2[pi:(pi+f) or None]]
                    pf = d.index(min(d)) #nearest forward point
                    pi += pf
    
                line.add(midpoint(p1[i],p2[pi]))
            if pi < len(p2)-1: #if last point of part 2 not included, extend
                line.add(midpoint(p1[i],p2[-1]))
            if (line[0].X==line[1].X) and (line[0].Y==line[1].Y):
                line.remove(0) #If first two points are identical, remove one
            #split line
    
            if line.count == 2:
                mid = midpoint(line[0],line[1])
                partA = arcpy.Array([line[0],mid])
                partB = arcpy.Array([mid,line[1]])
            elif line.count <2: print (line.count)
            else:
                partA = arcpy.Array([line[i] for i in range(line.count/2+line.count%2)])
                partB = arcpy.Array([line[i] for i in range(line.count/2,line.count)])
                
            iCur.insertRow([arcpy.Polyline(partA)])
            iCur.insertRow([arcpy.Polyline(partB)])
    
    del sCur, iCur

    #Tie centerlines to the nearest non-narrow segment of polygon
    arcpy.Near_analysis(centerline, extend_to, "", "LOCATION", "NO_ANGLE")
    
    uCur = arcpy.da.UpdateCursor(centerline,['SHAPE@','NEAR_X','NEAR_Y'])
    try:
        seg = uCur.next()
        while seg:
            partA = seg[0].getPart(0)
            if ((partA[0].X-seg[1])**2+(partA[0].Y-seg[2])**2)**.5 <= 18:
                partA.insert(0,arcpy.Point(seg[1],seg[2]))
            uCur.deleteRow()
            seg = uCur.next()
            partA.extend(seg[0].getPart(0))
            if ((partA[partA.count-1].X-seg[1])**2+ \
                (partA[partA.count-1].Y-seg[2])**2)**.5 <= 18:
                partA.add(arcpy.Point(seg[1],seg[2]))
            uCur.updateRow([arcpy.Polyline(partA),seg[1],seg[2]])
            seg = uCur.next()
    except StopIteration:
        return True

###******************* Part I *****************
#Step 0
arcpy.FeatureToLine_management(mupolygon, "MU_lines", "", "NO_ATTRIBUTES")
arcpy.Dissolve_management("MU_lines", "MU_lines_dis", "",\
                          "", "SINGLE_PART", "DISSOLVE_LINES")

arcpy.FeatureToPoint_management(mupolygon, "MU_point", "INSIDE")

arcpy.MakeFeatureLayer_management('MU_lines_dis', 'MU_lines_select')
arcpy.Delete_management('MU_lines_dis')
arcpy.Delete_management("MU_lines")
if query:
    arcpy.MakeFeatureLayer_management(mupolygon, 'MU_select', query)
    arcpy.SelectLayerByLocation_management("MU_lines_select","SHARE_A_LINE_SEGMENT_WITH",\
    'MU_select',"#","NEW_SELECTION")
    arcpy.SelectLayerByLocation_management("MU_lines_select","SHARE_A_LINE_SEGMENT_WITH",sapolygon,\
                                       '#',"REMOVE_FROM_SELECTION")
else:
    arcpy.SelectLayerByLocation_management("MU_lines_select","SHARE_A_LINE_SEGMENT_WITH",sapolygon,\
                                       '#',"NEW_SELECTION","INVERT")   

arcpy.SmoothLine_cartography("MU_lines_select", "MU_lines_gen", "PAEK",\
                             "25 Meters","FIXED_CLOSED_ENDPOINT", "NO_CHECK")
arcpy.Generalize_edit("MU_lines_gen", "2 Meters")
arcpy.SelectLayerByAttribute_management("MU_lines_select","SWITCH_SELECTION")
arcpy.Merge_management("MU_lines_select;MU_lines_gen", "MU_lines_gen_merge")

arcpy.FeatureToPolygon_management("MU_lines_gen_merge", "MU_gen",\
                                      "", "ATTRIBUTES", "MU_point")
arcpy.Delete_management("MU_lines_gen_merge")
arcpy.Delete_management("MU_point")
arcpy.Delete_management("MU_lines_gen")
if query:
    arcpy.MakeFeatureLayer_management('MU_gen', 'MU_gen_select', query)
    #Step 1
    #Expect to see several warnings of features dissappearing
    arcpy.Buffer_analysis("MU_gen_select", "MU_negbuff19", "-"+str(min_width)+" Meters")
        #Step 2
    arcpy.Buffer_analysis("MU_negbuff19", "MU_rebuff19", str(min_width+.1)+" Meters")
    
    #Step 3
    arcpy.Erase_analysis("MU_gen_select", "MU_rebuff19", \
                         "too_narrow", cluster_tolerance="0.1 Meters")
else:
    #Step 1
    #Expect to see several warnings of features dissappearing
    arcpy.Buffer_analysis("MU_gen", "MU_negbuff19", "-"+str(min_width)+" Meters")
    
    #Step 2
    arcpy.Buffer_analysis("MU_negbuff19", "MU_rebuff19",\
                          str(min_width+.1)+" Meters")
    
    #Step 3
    arcpy.Erase_analysis("MU_gen", "MU_rebuff19", "too_narrow", \
                         cluster_tolerance="0.1 Meters")

arcpy.Delete_management("MU_negbuff19")
arcpy.Delete_management("MU_rebuff19")

#Step 4
arcpy.MultipartToSinglepart_management("too_narrow", "too_narrow_sing")

#Step 5
arcpy.MakeFeatureLayer_management('too_narrow_sing','too_narrow_select',\
                                  'Shape_Length >=30 AND Shape_Length <=300')
    
#Step 6
arcpy.Buffer_analysis("too_narrow_select", "corner_buff_4", "0.4 Meters")

#Step 7 Intersect is necessary to produce a polygon county on both sides of narrows to equal 3
#Tolerance set to clean out polygons with 0 area
arcpy.Intersect_analysis("MU_gen #;corner_buff_4 #", \
                         "corner_inter", "ALL", "0.05 Meters")

#Step 8
arcpy.MultipartToSinglepart_management("corner_inter", "corner_inter_sing")

#Step 9
arcpy.SpatialJoin_analysis("corner_buff_4", "corner_inter_sing",\
            "corner_sum","JOIN_ONE_TO_ONE", "KEEP_ALL", "", "INTERSECT")

#Step 10
arcpy.MakeFeatureLayer_management('corner_sum','corner_sum_select',"Join_Count =2")

#Step 11
arcpy.SelectLayerByLocation_management('corner_sum_select',"CROSSED_BY_THE_OUTLINE_OF",\
                                       sapolygon,'#',"NEW_SELECTION",'INVERT')

#Step 12
arcpy.SelectLayerByLocation_management('too_narrow_select',"HAVE_THEIR_CENTER_IN"\
                                       ,'corner_sum_select','#',"NEW_SELECTION")

#Step 13
arcpy.Update_analysis("MU_gen", "too_narrow_select", "MU_cornered", "BORDERS", "0.1 Meters")

#Step 14
arcpy.MakeFeatureLayer_management("MU_cornered","MU_cornered_select")
arcpy.SelectLayerByLocation_management("MU_cornered_select","COMPLETELY_WITHIN",\
                                       'corner_sum_select')
arcpy.CopyFeatures_management("too_narrow_select",'truncated')
#Step 15 Eliminate
arcpy.Eliminate_management("MU_cornered_select", "MU_decornered", "LENGTH")

arcpy.Delete_management("too_narrow")
arcpy.Delete_management("too_narrow_sing")
arcpy.Delete_management("corner_buff_4")
arcpy.Delete_management("corner_inter")
arcpy.Delete_management("corner_inter_sing")
arcpy.Delete_management("corner_sum")
arcpy.Delete_management("MU_cornered")


###******************* Part II *****************


#*********************************************************
if query:
    arcpy.MakeFeatureLayer_management("MU_decornered", 'MU_decornered_select', query)
    #Step 1
    arcpy.Buffer_analysis("MU_decornered_select", "MU_2negbuff19","-19 Meters")
    
    #Step 2
    arcpy.Buffer_analysis("MU_2negbuff19", "MU_2rebuff19", "19.4 Meters")
    
    #Step 3
    arcpy.Erase_analysis('MU_decornered_select', "MU_2rebuff19", \
                         "too_narrow2", "0.1 Meters")
else:
    #Step 1
    arcpy.Buffer_analysis("MU_decornered", "MU_2negbuff19","-19 Meters")
    
    #Step 2
    arcpy.Buffer_analysis("MU_2negbuff19", "MU_2rebuff19", "19.4 Meters")
    
    #Step 3
    arcpy.Erase_analysis("MU_decornered", "MU_2rebuff19",\
                         "too_narrow2", "0.1 Meters")
arcpy.Delete_management("MU_2rebuff19")

#Step 4
arcpy.MultipartToSinglepart_management("too_narrow2", "too_narrow_sing2")

#Step 5
arcpy.FeatureToLine_management("too_narrow_sing2", "narrow_outlines")

#Step 6
arcpy.Intersect_analysis("narrow_outlines #;MU_decornered #", \
                         "narrow_outlines_inter", "ALL", "0.5 Meters")

#Step 7
arcpy.MakeFeatureLayer_management("narrow_outlines_inter","narrow_outlines_select",\
                                  "MUSYM <> MUSYM_1")
#Step 8 
arcpy.Dissolve_management("narrow_outlines_select", "casings", "FID_too_narrow_sing2")

#Step 9
arcpy.Densify_edit("casings", "DISTANCE", "30 Meters")

#Step 10 
centerline('casings','centerlines','MU_2negbuff19')

#Step 11
arcpy.Identity_analysis("centerlines", "MU_decornered", "centerlines_attributed", "ALL")

#Step 12
arcpy.Buffer_analysis("centerlines_attributed","inserts",\
                      "19 Meters", "FULL", "ROUND", "LIST", "FID_MU_decornered;MUSYM")

#Step 13
arcpy.MakeFeatureLayer_management('inserts','inserts_select')
arcpy.SelectLayerByLocation_management('inserts_select',"CROSSED_BY_THE_OUTLINE_OF",\
                                       sapolygon,'#',"NEW_SELECTION",'INVERT')
arcpy.CopyFeatures_management("inserts_select",'dilations')

#Step 14
arcpy.Update_analysis("MU_decornered", "dilations", "MU_inserted", \
                      "BORDERS", "0.1 Meters")

#Step 15
arcpy.Dissolve_management("MU_inserted", "MU_dilated", "MUSYM", "", \
                          "SINGLE_PART", "DISSOLVE_LINES")

arcpy.Delete_management("MU_2negbuff19")
arcpy.Delete_management("too_narrow2")
arcpy.Delete_management("too_narrow_sing2")
arcpy.Delete_management("narrow_outlines")
arcpy.Delete_management("narrow_outlines_inter")
arcpy.Delete_management("centerlines_attributed")
arcpy.Delete_management("inserts")
arcpy.Delete_management("MU_inserted")


### Densify 31.75 degrees to replace arcs with appropriate density
arcpy.Densify_edit("MU_dilated", "ANGLE", max_angle="31.75")

arcpy.FeatureToPoint_management("MU_decornered","MUpoints","INSIDE")

arcpy.SpatialJoin_analysis("MU_dilated", "MUpoints", "MU_final","JOIN_ONE_TO_ONE")
if query:
    arcpy.MakeFeatureLayer_management("MU_final", 'MU_final_select', query)
    arcpy.Buffer_analysis("MU_final_select", "MU_3negbuff19","-18.85 Meters")
    
    #Step
    arcpy.Buffer_analysis("MU_3negbuff19", "MU_3rebuff19", "19 Meters")
    
    #Step
    arcpy.Erase_analysis("MU_final_select", "MU_3rebuff19",\
                         "too_narrow3", "0.1 Meters")
else:
    arcpy.Buffer_analysis("MU_final", "MU_3negbuff19","-18.85 Meters")
    
    #Step
    arcpy.Buffer_analysis("MU_3negbuff19", "MU_3rebuff19", "19.25 Meters")
    
    #Step
    arcpy.Erase_analysis("MU_final", "MU_3rebuff19",\
                         "too_narrow3", "0.1 Meters")
    #Step
    arcpy.MultipartToSinglepart_management("too_narrow3", "Remaining_Narrows")

arcpy.Delete_management("MU_3negbuff19")
arcpy.Delete_management("MU_3rebuff19")
arcpy.Delete_management("too_narrow3")
arcpy.Delete_management("too_narrow3")

#https://www.ian-ko.com/ET_GeoWizards/WhitePapers/gw_smooth_generalize.htm
#http://desktop.arcgis.com/en/arcmap/10.3/analyze/modelbuilder/the-in-memory-workspace.htm#ESRI_SECTION1_1BDEBF62493E489FA7A1CD7E4E951D5A