import arcpy import math def JZpointmake(inpath, outfileshp): # region [读取shp] shpstr = 'JZD1.shp' outpath = outfileshp + '\\' + 'JZD1.shp' try: arcpy.CreateFeatureclass_management(outfileshp, shpstr, 'POINT') except: arcpy.Delete_management(outpath) arcpy.CreateFeatureclass_management(outfileshp, shpstr, 'POINT') # 加载 arcpy.MakeFeatureLayer_management(inpath, 'LB') arcpy.MakeFeatureLayer_management(outpath, 'JZD') sr = arcpy.Describe('LB').spatialReference # 添加一群字段 arcpy.AddField_management('JZD', '标识码', 'LONG') arcpy.AddField_management('JZD', 'YSDM', 'TEXT') arcpy.AddField_management('JZD', 'JZDH', 'TEXT') arcpy.AddField_management('JZD', 'JBLX', 'TEXT') arcpy.AddField_management('JZD', 'JZDLX', 'TEXT') arcpy.AddField_management('JZD', 'ZDBH', 'TEXT') arcpy.AddField_management('JZD', '顺序号', 'LONG') arcpy.AddField_management('JZD', 'flag', 'LONG') #float按理说应该符合的 bsm = 0 field1 = ['SHAPE@', '标识码', 'YSDM', 'JZDH', 'JBLX', 'JZDLX', 'ZDBH', '顺序号','flag'] # endregion with arcpy.da.InsertCursor('JZD', field1) as cursorZ: with arcpy.da.SearchCursor('LB', ['SHAPE@', '宗地编号']) as cursorA: for curA in cursorA: # region [读折点] poly1 = curA[0] zdbh = curA[1] if zdbh == None or zdbh == '' or zdbh == ' ': pass else: # 获取折点 points = poly1.getPart() pt = points[0] mm = 0 linelist = [] indexp = 0 #存的都是点 while mm < pt.count: jzd = pt[mm] try: xx = round(jzd.X,4) yy = round(jzd.Y,4) if mm == 0 or mm == indexp: xx0 = xx yy0 = yy slist = [] slist.append(pt[mm]) else: numdis0 = xx - xx0 numdis1 = yy - yy0 numdis = math.sqrt((numdis0 * numdis0)+(numdis1*numdis1)) if xx == xx0 and yy == yy0: #一个的结束 indexp = mm + 1 linelist.append(slist) if numdis < 0.1: pass else: slist.append(pt[mm]) mm = mm + 1 except: if mm == indexp: indexp = mm + 1 mm = mm + 1 llen = len(linelist) #判断下如果多环，排个序 linelist1 = pdlen(linelist) linelist = linelist1 # endregion #甭管有没有环，都是第一个来处理 ln1 = linelist[0] list1 = [] nn1 = 0 while nn1 < len(ln1): if nn1 == len(ln1) - 1: fpx = round(ln1[nn1].X,4) fpy = round(ln1[nn1].Y,4) lpx = round(ln1[0].X,4) lpy = round(ln1[0].Y,4) fpoint = arcpy.Point(fpx, fpy) lpoint = arcpy.Point(lpx, lpy) else: fpx = round(ln1[nn1].X,4) fpy = round(ln1[nn1].Y,4) lpx = round(ln1[nn1 + 1].X,4) lpy = round(ln1[nn1 + 1].Y,4) fpoint = arcpy.Point(fpx, fpy) lpoint = arcpy.Point(lpx, lpy) new_array = arcpy.Array() new_array.append(fpoint) new_array.append(lpoint) ln2 = arcpy.Polyline(new_array,sr) #先把起始点加进去 bl1 = addpoint(list1,fpoint) if bl1 == 1: list1.append(fpoint) # list1.append(lpoint) #插入点 arcpy.SelectLayerByLocation_management('LB','INTERSECT',ln2) count0 = int(arcpy.GetCount_management('LB').getOutput(0)) if count0 < 2: bl4 = addpoint(list1,lpoint) if bl4 == 1: list1.append(lpoint) nn1 = nn1 + 1 else: #对于同一个线段来说，把所有点的位置准确点 listp1 = [] with arcpy.da.SearchCursor('LB', ['SHAPE@', '宗地编号']) as cursorB: for curB in cursorB: zdbh1 = curB[1] #排除自己 if zdbh1 == zdbh: pass else: try: shpB = curB[0] lineB = shpB.boundary () rln = lineB.intersect (ln2, 2) fpoint1 = rln.firstPoint lpoint1 = rln.lastPoint fx = round(fpoint1.X,4) fy = round(fpoint1.Y,4) lx = round(lpoint1.X,4) ly = round(lpoint1.Y,4) #判断这两个点是否都在里面 fp = arcpy.Point(fx, fy) lp = arcpy.Point(lx, ly) list2 = [] bl1 = addpoint(list1,fp) if bl1 == 1: #判断与起点的距离 dis1 = math.sqrt(((fpx - fx)*(fpx - fx)) + ((fpy-fy)*(fpy-fy))) slist2 = [] slist2.append(dis1) slist2.append(fp) list2.append(slist2) bl1 = addpoint(list1,lp) if bl1 == 1: dis1 = math.sqrt(((fpx - lx)*(fpx - lx)) + ((fpy-ly)*(fpy-ly))) slist2 = [] slist2.append(dis1) slist2.append(lp) list2.append(slist2) if len(list2) == 1: listp1.append(list2[0][1]) elif len(list2) == 2: #对2排序 if list2[0][0] < list2[1][0]: listp1.append(list2[0][1]) listp1.append(list2[1][1]) else: listp1.append(list2[1][1]) listp1.append(list2[0][1]) except: pass #得到所有的点，进行一个与起点的排序 nn2 = 0 listp2 =[] while nn2 < len(listp1): num0 = fpx - listp1[nn2].X num1 = fpy - listp1[nn2].Y dis2 = math.sqrt((num0*num0)+(num1*num1)) slistp2 = [] slistp2.append(dis2) slistp2.append(listp1[nn2]) listp2.append(slistp2) nn2 = nn2 + 1 listp2.sort() nn3 = 0 bl3 = 0 while nn3 < len(listp2): #这里要看下要不要除重 bl3 = addpoint(list1,listp2[nn3][1]) if bl3 == 1: list1.append(listp2[nn3][1]) nn3 = nn3 + 1 bl4 = addpoint(list1,lpoint) if bl4 == 1: list1.append(lpoint) nn1 = nn1 + 1 if len(linelist) == 1: listF = linezd(list1,poly1) #直接写 ii = 0 zdh = 0 while ii < len(listF): if ii == 0: zdh = zdh + 1 bsm = bsm + 1 # shape point1 = arcpy.Point(listF[ii][0], listF[ii][1]) # 要素代码 ysdm = '6002010300' # 顺序号 sxh = ii + 1 # 界址点号 jzdh = 'J' + str(zdh) # 宗地代码 zddm = curA[1] rlist = [] rlist.append(point1) rlist.append(bsm) rlist.append(ysdm) rlist.append(jzdh) rlist.append('') rlist.append('') rlist.append(zddm) rlist.append(zdh) rlist.append(1) cursorZ.insertRow(rlist) else: numdis0 = listF[ii][0] - listF[ii - 1][0] numdis1 = listF[ii][1] - listF[ii - 1][1] numdis2 = math.sqrt((numdis0 * numdis0)+(numdis1*numdis1)) if numdis2 < 0.1: pass else: zdh = zdh + 1 bsm = bsm + 1 # shape point1 = arcpy.Point(listF[ii][0], listF[ii][1]) # 要素代码 ysdm = '6002010300' # 顺序号 sxh = ii + 1 # 界址点号 jzdh = 'J' + str(zdh) # 宗地代码 zddm = curA[1] rlist = [] rlist.append(point1) rlist.append(bsm) rlist.append(ysdm) rlist.append(jzdh) rlist.append('') rlist.append('') rlist.append(zddm) rlist.append(zdh) rlist.append(1) cursorZ.insertRow(rlist) ii = ii + 1 else: mm0 = 0 zdh = 0 while mm0 < len(linelist): if mm0 == 0: listF = linezd(list1,poly1) ii = 0 while ii < len(listF): if ii == 0: zdh = zdh + 1 bsm = bsm + 1 # shape point1 = arcpy.Point(listF[ii][0], listF[ii][1]) # 要素代码 ysdm = '6002010300' # 顺序号 sxh = ii + 1 # 界址点号 jzdh = 'J' + str(zdh) # 宗地代码 zddm = curA[1] rlist = [] rlist.append(point1) rlist.append(bsm) rlist.append(ysdm) rlist.append(jzdh) rlist.append('') rlist.append('') rlist.append(zddm) rlist.append(zdh) rlist.append(mm0+1) cursorZ.insertRow(rlist) else: numdis0 = listF[ii][0] - listF[ii - 1][0] numdis1 = listF[ii][1] - listF[ii - 1][1] numdis2 = math.sqrt((numdis0 * numdis0)+(numdis1*numdis1)) if numdis2 < 0.1: pass else: zdh = zdh + 1 bsm = bsm + 1 # shape point1 = arcpy.Point(listF[ii][0], listF[ii][1]) # 要素代码 ysdm = '6002010300' # 顺序号 sxh = ii + 1 # 界址点号 jzdh = 'J' + str(zdh) # 宗地代码 zddm = curA[1] rlist = [] rlist.append(point1) rlist.append(bsm) rlist.append(ysdm) rlist.append(jzdh) rlist.append('') rlist.append('') rlist.append(zddm) rlist.append(zdh) rlist.append(mm0+1) cursorZ.insertRow(rlist) ii = ii + 1 else: listF = linezd(linelist[mm0],poly1) ii = len(listF) - 1 while ii >= 0: if ii == len(listF) - 1: zdh = zdh + 1 bsm = bsm + 1 # shape point1 = arcpy.Point(listF[ii][0], listF[ii][1]) # 要素代码 ysdm = '6002010300' # 顺序号 sxh = ii + 1 # 界址点号 jzdh = 'J' + str(zdh) # 宗地代码 zddm = curA[1] rlist = [] rlist.append(point1) rlist.append(bsm) rlist.append(ysdm) rlist.append(jzdh) rlist.append('') rlist.append('') rlist.append(zddm) rlist.append(zdh) rlist.append(mm0+1) cursorZ.insertRow(rlist) else: numdis0 = listF[ii][0] - listF[ii - 1][0] numdis1 = listF[ii][1] - listF[ii - 1][1] numdis2 = math.sqrt((numdis0 * numdis0)+(numdis1*numdis1)) if numdis2 < 0.1: pass else: zdh = zdh + 1 bsm = bsm + 1 # shape point1 = arcpy.Point(listF[ii][0], listF[ii][1]) # 要素代码 ysdm = '6002010300' # 顺序号 sxh = ii + 1 # 界址点号 jzdh = 'J' + str(zdh) # 宗地代码 zddm = curA[1] rlist = [] rlist.append(point1) rlist.append(bsm) rlist.append(ysdm) rlist.append(jzdh) rlist.append('') rlist.append('') rlist.append(zddm) rlist.append(zdh) rlist.append(mm0+1) cursorZ.insertRow(rlist) ii = ii - 1 mm0 = mm0 + 1 arcpy.AddMessage(zdbh + ' FINISH') def pdlen(list1): ii = 0 list2 = [] while ii < len(list1): pts = list1[ii] length1 = 0 mm = 0 while mm < len(pts): if mm == len(pts) - 1: pt1 = pts[mm] pt2 = pts[0] dis1 = math.sqrt(((pt1.X-pt2.X)*(pt1.X-pt2.X))+((pt1.Y-pt2.Y)*(pt1.Y-pt2.Y))) length1 = length1 + dis1 mm = mm + 1 else: pt1 = pts[mm] pt2 = pts[mm + 1] dis1 = math.sqrt(((pt1.X-pt2.X)*(pt1.X-pt2.X))+((pt1.Y-pt2.Y)*(pt1.Y-pt2.Y))) length1 = length1 + dis1 mm = mm + 1 slist2 = [] slist2.append(length1) slist2.append(pts) list2.append(slist2) ii = ii + 1 list2.sort(reverse=True) nn = 0 list3 = [] while nn < len(list2): list3.append(list2[nn][1]) nn = nn + 1 return list3 def addpoint(list1,pt): ptx = round(pt.X,4) pty = round(pt.Y,4) ii = 0 bl = 1 while ii < len(list1): if ptx == round(list1[ii].X,4) and pty == round(list1[ii].Y,4): bl = -1 ii = ii + 1 return bl def linezd(line1,poly): # 获取折点 mm = 0 pp = 0 list1 = [] while mm < len(line1): plist = [] jzd = line1[mm] try: xx = jzd.X yy = jzd.Y plist.append(xx) plist.append(yy) # 初始顺序 plist.append(pp) list1.append(plist) pp = pp + 1 mm = mm + 1 except: mm = mm + 1 # 对折点排序（算法） if len(list1) < 3: listF = [] else: # listF = zdsort(list1, zx) listF = jzdsort(list1,poly) return listF def jzdsort(list0,poly): m = 0 n = 0 northwest = [] northwest.append(poly.extent.XMin) northwest.append(poly.extent.YMax) minDis = (northwest[1] - list0[0][1])*(northwest[1] - list0[0][1]) + (northwest[0] - list0[0][0])*(northwest[0] - list0[0][0]) maxStdangle = cal_stdangle(northwest,list0[0]) while m < len(list0): Dis = (northwest[1] - list0[m][1])*(northwest[1] - list0[m][1]) + (northwest[0] - list0[m][0])*(northwest[0] - list0[m][0]) if Dis < minDis: minDis = Dis maxStdangle = cal_stdangle(northwest,list0[m]) n = m if abs(Dis-minDis) < 0.000001: stdangle = cal_stdangle(northwest,list0[m]) if stdangle > maxStdangle: maxStdangle = stdangle n = m m = m + 1 i = 0 while i < len(list0): if n == 0: vectangleN = cal_stdangle(list0[0], list0[1]) vectangleNl = cal_stdangle(list0[len(list0)-1], list0[0]) if abs(vectangleN-vectangleNl) < 3: n = len(list0)-1 else: if n == len(list0)-1: vectangleN = cal_stdangle(list0[n], list0[0]) else: vectangleN = cal_stdangle(list0[n], list0[n+1]) vectangleNl = cal_stdangle(list0[n-1], list0[n]) if abs(vectangleN-vectangleNl) < 3: n = n-1 i = i + 1 list4 = [] a = n num = 0 while a < len(list0): slist4 = [] num = num + 1 slist4.append(list0[a][0]) slist4.append(list0[a][1]) slist4.append(num) list4.append(slist4) a = a + 1 a = 0 while a < n: slist4 = [] num = num + 1 slist4.append(list0[a][0]) slist4.append(list0[a][1]) slist4.append(num) list4.append(slist4) a = a + 1 return list4 def cal_stdangle(lastp,p): stdangle=0 dy= p[1]-lastp[1] dx= p[0]-lastp[0] if dx==0 and dy>0: stdangle = 0 if dx==0 and dy<0: stdangle = 180 if dy==0 and dx>0: stdangle = 90 if dy==0 and dx<0: stdangle = 270 if dx>0 and dy>0: stdangle = math.atan(dx/dy)*180/math.pi elif dx<0 and dy>0: stdangle = 360 + math.atan(dx/dy)*180/math.pi elif dx<0 and dy<0: stdangle = 180 + math.atan(dx/dy)*180/math.pi elif dx>0 and dy<0: stdangle = 180 + math.atan(dx/dy)*180/math.pi return stdangle def calcangle(lastp,p): angle=0 stdangle = cal_stdangle(lastp,p) if (stdangle < 45 and stdangle >= 0) or (stdangle <= 360 and stdangle >= 315): angle = 270 if (stdangle < 135 and stdangle >= 45): angle = 0 if (stdangle < 225 and stdangle >= 135): angle = 90 if (stdangle < 315 and stdangle >= 225): angle = 180 return angle if __name__ == '__main__': try: # 输入shp，输出界址点文件夹 JZpointmake(r'D:\2work_now\20221026LQ\1231\Export_Output.shp',r'D:\2work_now\20221026LQ\1231\JZD') except arcpy.ExecuteError: arcpy.AddMessage(arcpy.GetMessages())