面的节点夹角大于30并且小于150,距离外界矩形左上角最近的点为第一个点,对于多部件每一个都修改,修改面的开始点位置,运行界面如图14-10所示。
#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
标签:return,pt,Python,ArcGIS,y1,左上角,y2,arcpy,pt2 From: https://www.cnblogs.com/suoyike1001/p/16998584.html