導航

#coding=utf8
import arcpy

import os
import sys
import math
isEdit=False
def getdis(x1,y1,x2,y2): #獲得兩個點距離
    return math.sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))
###############
def getXYAngle(x1,y1,x2,y2):#點x1,y1和點x2,y2的角度
    if (math.fabs(x1 - x2) < 0.00001):
        if (y1 < y2):
            return 90
        else:
            return 270;
    else:
        k = (y2 - y1) / (x2 - x1);
        k = math.atan(k) * 180 / math.pi
        if (x2 < x1):  #二、三象限
            return k + 180;
        elif (y2 >= y1): #//一象限
            return k;
        else:#四象限
            return k + 360

def getAngle(pt1,pt2):#獲得兩個點的角度
    x1=pt1.X
    y1=pt1.Y
    x2=pt2.X
    y2=pt2.Y
    return getXYAngle(x1,y1,x2,y2)



def getLineAngle(pt1, pt2, pt3,mydir):#mydir是true順時針,false是逆時針
    if (pt1==None):
        #arcpy.AddMessage("pt1為空")
        return -1
    if (pt2==None):
        #arcpy.AddMessage("pt2為空")
        return -2
    if (pt3==None):
        #arcpy.AddMessage("pt3為空")
        return  -3

    A1 = getAngle(pt1, pt2)
    A2 = getAngle(pt2, pt3)
    angle = 0;
    if mydir:
        angle = 180 + A2 - A1
    else:
        angle = 180 + A1 - A2
    if (angle >= 360):
        angle = angle - 360
    if (angle < 0):
        angle = angle + 360
    return angle


def getintersectionAnge(pt1, pt2, pt3):#獲得夾角
    a=getLineAngle(pt1, pt2, pt3,True)
    if a>180:
        a=a-180
    return a
def getminXmaxY(array):#獲得最小的X和最大Y
    num=len(array)
    minx=999999999999
    maxy=0;

    for i in range(num):
        pt=array[i]
        if pt:
            if pt.X<minx:
                minx=pt.X
            </minx:
if pt.Y>maxy:
                maxy=pt.Y
    return minx,maxy

###########
def splitNgeometry(mgeometry):
    num=mgeometry.count
    Sumarray = arcpy.Array()
    parray = arcpy.Array()
    for i in range(num): #pntcount
        pt=mgeometry[i]
        if pt:
            parray.add(pt)
        else:#內邊形
            Sumarray.add(parray)
            parray.removeAll()
    Sumarray.add(parray)
    return Sumarray

################
def MINPoint(partgeometry):#按照左上點，修改圖形
    global isEdit

    Topx,Topy=getminXmaxY(partgeometry)
    num=partgeometry.count

    #arcpy.AddMessage("partgeometry.extent："+str(Topx)+":y="+str(Topy))
    maxd=99999999999;
    idx=0;
    for i in range(num):

        pt=partgeometry[i]

        x=pt.X
        y=pt.Y

        d= getdis(x,y,Topx,Topy)
        #arcpy.AddMessage("KKK"+str(i)+",坐標x="+str(x)+",y="+str(y)+",d="+str(d))
        if (d<maxd):
            </maxd):
#計算夾角
            if (i>1):
                p1=partgeometry[i-1]
            else:
                p1=partgeometry[num-1]

            if i<(num-1):
                p2=partgeometry[i+1]
            else:
                p2=partgeometry[0]
            myAngle= getintersectionAnge(p1, pt, p2)

            #arcpy.AddMessage("i:"+str(i)+"，myAngle="+str(myAngle))

            if myAngle>=30 and myAngle<=150:
                idx=i
                maxd=d
    #arcpy.AddMessage("idx============:"+str(idx)+"，min="+str(maxd))
    if idx<1:

        return partgeometry
    else:
        isEdit=True

        array = arcpy.Array()
        for i in range(idx,num):
            #arcpy.AddMessage("i===========:"+str(i))
            if partgeometry[i]:

                array.add(partgeometry[i])
        for i in range(idx):
            #arcpy.AddMessage("i++++++++++:"+str(i))
            if partgeometry[i]:
                array.add(partgeometry[i])
        return array

inFeature  = arcpy.GetParameterAsText(0)
outFeature  = arcpy.GetParameterAsText(1)
arcpy.Select_analysis(inFeature, outFeature, '')
desc = arcpy.Describe(outFeature)
shapeName = desc.ShapeFieldName
OIDField=desc.OIDFieldName

result = arcpy.GetCount_management(inFeature)
count= int(result.getOutput(0))
if count < 1:
    arcpy.AddMessage(inFeature+u"沒有數據")
else:

    rows = arcpy.UpdateCursor(outFeature)

    n=1

    try:
        for row in rows:
            isEdit=False
            arcpy.SetProgressorPosition()
            arcpy.SetProgressorLabel(u"正在等待，完成"+str(round(1.0*n*100/count,1))+"%...")
            geometry = row.getValue(shapeName)
            FID=row.getValue(OIDField)
            part_count = geometry.partCount #有幾部分
            #arcpy.AddMessage("FID:"+str(FID)+",part_count:"+str(part_count))
            Sumarray = arcpy.Array()
            for i in range(part_count):
                partgeometry=geometry.getPart(i)
                SpliArray=splitNgeometry(partgeometry)
                N=SpliArray.count
                #arcpy.AddMessage("NNNNN=====:"+str(N))
                for j in range(N):
                    Splitgeometry=SpliArray[j]

                    array=MINPoint(Splitgeometry)
                #if (part_count>1):
                    try:
                        Sumarray.add(array)
                    except Exception as err:
                        arcpy.AddError(u"錯誤=============j:"+str(j)+","+err.message)

            if isEdit:
                #row.setValue(shapeName, array)
                row.setValue(shapeName, Sumarray)
                rows.updateRow(row)
                arcpy.AddMessage("FID:"+str(FID)+u"修改")
            n=n+1

    finally:
        arcpy.ResetProgressor()
        #del row
        del rows