14.1 表属性中一个字出现次数计算
#coding=utf8
import arcpy from arcpy import env import os import sys import time import ylpy import numpy def getidx(List,mystr):#获得列表的中位置 try: idx=List.index(mystr) return idx except Exception as e: # 都能处理异常 return -1 fc= arcpy.GetParameterAsText(0) fieldname= arcpy.GetParameterAsText(1) tablename= arcpy.GetParameterAsText(2) num=ylpy.getCount(fc) ylpy.initProgress("get",num) rows = arcpy.da.SearchCursor(fc,(fieldname)) start = time.clock() nameList=[] nameint=[] ######################################### ## for row in rows: name =u""+row[0].strip() for s in name: mystr=s #decode('utf-8') idx=getidx(nameList,mystr) if (idx>-1): nameint[idx]=nameint[idx]+1 else: arcpy.AddMessage(u"______"+mystr) nameint.append(1) nameList.append(mystr) ylpy.step() ylpy.freeProgress() del rows mynumpy=numpy.rec.fromarrays([nameList,nameint],names=["name","n"]) arcpy.da.NumPyArrayToTable(mynumpy, tablename) elapsed = (time.clock() - start) arcpy.AddMessage("Time used:"+str(elapsed)) #获得运行时间
14.2 矢量数据批量裁剪
#coding=utf8 import sys, os, string import arcpy from arcpy import env def getuniqueValue(inTable,inField):#获得字段唯一值 rows = arcpy.SearchCursor(inTable) # Create an empty list uniqueList = [] try: for row in rows: # If the value is not already in the list, append it if row.getValue(inField) not in uniqueList: uniqueList.append(row.getValue(inField)) return uniqueList finally: if row: del row if rows: del rows def getlength(mystr):#获得长度，汉字长度为2，英文为1 lenTxt = len(mystr) lenTxt_utf8 = len(mystr.decode('utf-8')) size = int((lenTxt_utf8 - lenTxt)/2 + lenTxt) #arcpy.AddMessage("'{0}'最后长度{1},lenTxt={2},lenTxt_utf8={3}".format(mystr,size,lenTxt,lenTxt_utf8)) return size def midFill(sumn,mystr,Fill):#填充长度一样 n=getlength(mystr) if n>=sumn: return mystr leftn=int((sumn-n)/2) s="" lefts=s.ljust(leftn,Fill) s="" rightn=sumn-n-leftn rights=s.ljust(rightn,Fill) return lefts+mystr+rights def printauthor(toolname):#打印作者信息 titlestr="" sumn=60 Fill='*' titlestr=titlestr.ljust(sumn,Fill) arcpy.AddMessage(titlestr) arcpy.AddMessage(midFill(sumn,u"欢迎使用:"+toolname,Fill)) mystr=u"本工具闫磊编制,QQ:276529800" arcpy.AddMessage(midFill(sumn,mystr,Fill)) mystr=u"使用前请做好数据备份，工具产生的不良后果请自行承担！" arcpy.AddMessage(midFill(sumn,mystr,Fill)) arcpy.AddMessage(titlestr) arcpy.env.overwriteOutput = True inworkspace = arcpy.GetParameterAsText(0) clipshp = arcpy.GetParameterAsText(1) fieldname= arcpy.GetParameterAsText(2) outworkspace = arcpy.GetParameterAsText(3) gdbbool = arcpy.GetParameterAsText(4).upper() desc = arcpy.Describe(clipshp) shapeName = desc.shapeFieldName #shape字段 result = arcpy.GetCount_management(clipshp) count= int(result.getOutput(0)) if count <= 0: arcpy.AddMessage(clipshp+u"没有数据") pass printauthor(u"矢量数据批量裁剪") uniqueList=getuniqueValue(clipshp,fieldname) Dissolveb=False;#是否融合 num1=len(uniqueList) if not count==num1: arcpy.AddMessage(u"由于"+fieldname+u"字段不是唯一值,软件做了融合处理，你看到几个和最终结果个数不一致,原始有"+str(count)+u"个,最后输出只有"+str(num1)+u"个数据库") outnewshp = arcpy.CreateUniqueName("yl_temp") #临时 arcpy.Dissolve_management(clipshp, outnewshp, [fieldname], "", "MULTI_PART","DISSOLVE_LINES") count=len(uniqueList) clipshp= outnewshp Dissolveb=True mydatasets= inworkspace.split(";") num=len(mydatasets)*count arcpy.SetProgressor("step", u"更新"+fieldname,0,num,1) rows = arcpy.SearchCursor(clipshp) row = rows.next() i=0; while row: #arcpy.AddMessage(u"7=执行到这里") fieldvalue =""+ str(row.getValue(fieldname)) out_mdb="" #arcpy.AddMessage("======================================================out_mdb"+out_mdb) if gdbbool=="MDB": out_mdb=outworkspace + "\"+fieldvalue+".mdb" #os.path.basename(dataset) elif gdbbool=="GDB": out_mdb=outworkspace + "\"+fieldvalue+".gdb" else:#shp out_mdb=outworkspace +os.sep+fieldvalue arcpy.AddMessage(u"out_mdb"+out_mdb) if not arcpy.Exists(out_mdb): #可以文件夹 if gdbbool=="MDB": arcpy.CreatePersonalGDB_management(os.path.dirname(out_mdb),os.path.basename(out_mdb)) elif gdbbool=="GDB": arcpy.CreateFileGDB_management(os.path.dirname(out_mdb),os.path.basename(out_mdb)) else:#shp os.makedirs(out_mdb) #arcpy.AddMessage("88888888888888888888888888888888888888") geometry=row.getValue(shapeName) for dataset in mydatasets: i+=1 arcpy.SetProgressorPosition() arcpy.SetProgressorLabel(u"正在裁剪....,完成："+str(i*100/num)+"%" ) try: #mylayer=os.path.basename(dataset) #原来有别名的裁剪有问题 desc=arcpy.Describe(dataset) mylayer=desc.baseName arcpy.AddMessage(u"clip:"+dataset+" to "+out_mdb+"\"+ mylayer) mylayer=mylayer.replace("(","") mylayer=mylayer.replace(")","") #oldoutFC=out_mdb+"\"+ mylayer outFC = arcpy.ValidateTableName(mylayer,out_mdb) #arcpy.AddMessage("outFC:"+outFC) #arcpy.Clip_analysis(dataset, jfb_Select,out_mdb+"\"+ mylayer, "") if gdbbool=="SHP": arcpy.Clip_analysis(dataset, geometry,out_mdb+"\"+outFC+".shp", "") else: arcpy.Clip_analysis(dataset, geometry,out_mdb+"\"+outFC, "") except Exception, ErrorDesc: #If an error set output boolean parameter "Error" to True. arcpy.AddError(str(ErrorDesc)) row = rows.next() if Dissolveb: if arcpy.Exists(clipshp): arcpy.Delete_management(clipshp) if row: del row arcpy.ResetProgressor()
14.3 矢量数据批量合库
#coding=utf8 import sys ############# ################################# import arcpy import string import os def isGDB(yldir): if (yldir.lower().endswith(".gdb")): return True elif (yldir.lower().endswith(".mdb")): return True else: return False def CopyDir(outdb,workspace): #shp arcpy.env.workspace = workspace files = arcpy.ListWorkspaces("","") if files: for File in files: if (File==outdb): continue AppendGDB(File,outdb) def AppendGDB(File,outdb): if (File.lower()==outdb.lower()): return arcpy.AddMessage('File=========:'+File) arcpy.env.workspace = outdb fcs = arcpy.ListFeatureClasses() for fc in fcs: try: if arcpy.Exists(File + "\" + fc): arcpy.AddMessage("fc:"+File + "\" + fc) arcpy.Append_management([ File + "\" + fc], outdb + "\" + fc,"NO_TEST","","") elif arcpy.Exists(File + os.sep + fc + ".shp"): # 或者shp: arcpy.AddMessage("fc:" + File + "\" + fc+".shp") arcpy.Append_management([File + "\" + fc+".shp"], outdb + "\" + fc, "NO_TEST", "", "") else: arcpy.AddMessage("not exists:"+File + "\" + fc) except arcpy.ExecuteError: arcpy.AddWarning(arcpy.GetMessages()) fcs = arcpy.ListTables() for fc in fcs: arcpy.AddMessage("fc:"+fc) try: if arcpy.Exists(File + os.sep + fc):#或者dbf arcpy.Append_management([File + "\" + fc], outdb + "\" + fc,"NO_TEST","","") elif arcpy.Exists(File + os.sep + fc+".dbf"): arcpy.Append_management([File + "\" + fc+".dbf"], outdb + "\" + fc, "NO_TEST", "", "") else: arcpy.AddMessage("not exists:"+File + "\" + fc) except arcpy.ExecuteError: arcpy.AddWarning(arcpy.GetMessages()) dss = arcpy.ListDatasets() for ds in dss: arcpy.AddMessage("ds:"+ds) arcpy.env.workspace = outdb+"\"+ds fcs1 = arcpy.ListFeatureClasses() for fc1 in fcs1: arcpy.AddMessage("fc1:"+fc1) try: if arcpy.Exists(File + "\" + ds + "\" + fc1): arcpy.Append_management([File + "\" + ds + "\" + fc1], outdb + "\" + ds + "\" + fc1,"NO_TEST","","") else: arcpy.AddMessage("not exists:"+File + "\" + ds + "\" + fc1) except arcpy.ExecuteError: arcpy.AddWarning(arcpy.GetMessages()) workspace =arcpy.GetParameterAsText(0) #'C:UsersAdministratorDesktop\cc' outdb =arcpy.GetParameterAsText(1) #'C:UsersAdministratorDesktop\lutian.mdb' for dirpath, dirnames, filenames in os.walk(workspace): arcpy.AddMessage('dirpath=======:'+dirpath) if isGDB(dirpath):#是gdb arcpy.AddMessage(u'dirpath=======是gdb:'+dirpath) AppendGDB(dirpath,outdb) for dirname in dirnames: arcpy.AddMessage("dirname=="+dirname) if isGDB(dirname): arcpy.AddMessage(u'dirname=======是gdb:'+dirpath) continue else: filepath= os.path.join(dirpath,dirname) AppendGDB(filepath,outdb) arcpy.AddMessage(u'dirname不是gdb:'+dirname) for filename in filenames: if filename.lower().endswith(".mdb"): filepath= os.path.join(dirpath,filename) arcpy.AddMessage('filepath===:'+filepath) AppendGDB(filepath,outdb)
14.4 影像批量裁剪
#coding=utf8 import sys, os, string,types import arcpy from arcpy import env def getuniqueValue(inTable,inField): rows = arcpy.SearchCursor(inTable) # Create an empty list uniqueList = [] try: for row in rows: # If the value is not already in the list, append it if row.getValue(inField) not in uniqueList: uniqueList.append(row.getValue(inField)) return uniqueList finally: if row: del row if rows: del rows arcpy.env.overwriteOutput = True oldraster = arcpy.GetParameterAsText(0) clipshp = arcpy.GetParameterAsText(1) fieldname= arcpy.GetParameterAsText(2) outworkspace= arcpy.GetParameterAsText(3) outdesc = arcpy.Describe(outworkspace) ext= arcpy.GetParameterAsText(4) if outdesc.dataType== "Workspace": ext="" elif not outdesc.dataType=="Folder":#如FeatureDataset FeatureLayer arcpy.AddError(u"格式错误无法裁剪") pass arcpy.CheckOutExtension("spatial") desc = arcpy.Describe(clipshp) shapeName = desc.shapeFieldName #shape字段 result = arcpy.GetCount_management(clipshp) num= int(result.getOutput(0)) if num <= 0: arcpy.AddMessage(clipshp+u"没有数据") pass uniqueList=getuniqueValue(clipshp,fieldname) Dissolveb=False;#是否融合 num1=len(uniqueList) if not num==num1: arcpy.AddMessage(u"由于"+fieldname+u"字段不是唯一值,软件做了融合处理，" +"你看到几个和最终结果个数不一致,原始有" +str(num)+u"个,最后输出只有"+str(num1)+u"个") outnewshp = arcpy.CreateUniqueName("yl_temp") #临时 arcpy.Dissolve_management(clipshp, outnewshp, [fieldname], "", "MULTI_PART","DISSOLVE_LINES") num=len(uniqueList) clipshp= outnewshp Dissolveb=True arcpy.SetProgressor("step", u"正在裁剪",0,num,1) rows = arcpy.SearchCursor(clipshp) i=0 for row in rows: try: i+=1 arcpy.SetProgressorPosition() arcpy.SetProgressorLabel(u"正在裁剪....,完成："+str(i*100/num)+"%" ) fieldvalue=str(row.getValue(fieldname)) if fieldvalue==None: fieldvalue="None" geometry=row.getValue(shapeName) if (outdesc.dataType== "Workspace"): outFC = arcpy.ValidateTableName(fieldvalue+ext,outworkspace) out_raster =outworkspace+"/"+outFC else: out_raster=outworkspace+"/"+fieldvalue+ext arcpy.AddMessage(u"正在裁剪:"+out_raster); #out_raster =outworkspace+"/"+fieldvalue+ext #arcpy.sa.ExtractByMask(oldraster, geometry, out_raster) arcpy.Clip_management(oldraster,"#",out_raster,geometry, "#", "ClippingGeometry") arcpy.env.pyramid = "PYRAMIDS 3 BILINEAR JPEG" arcpy.BuildPyramids_management(out_raster) except Exception, ErrorDesc: #If an error set output boolean parameter "Error" to True. arcpy.AddError(str(ErrorDesc)) arcpy.ResetProgressor() if Dissolveb: if arcpy.Exists(clipshp): arcpy.Delete_management(clipshp) if row: del row if rows: del rows
11.5 按数据库标准创建要素类和表
#coding=utf8 import arcpy import os import sys import math from arcpy.sa import * #初始化进度条 def initProgress(hint,num): arcpy.SetProgressor("step", u"正在"+hint,0,num,1) #进度条 def step(): arcpy.SetProgressorLabel(u"正在进行....") arcpy.SetProgressorPosition() #释放进度条 def freeProgress(): arcpy.ResetProgressor() def getCount(inFeature): result = arcpy.GetCount_management(inFeature) count= int(result.getOutput(0)) return count #获得唯一值 def getuniqueValue(inTable,inField): rows = arcpy.da.SearchCursor(inTable,inField) # Create an empty list uniqueList = [] try: for row in rows: # If the value is not already in the list, append it if row[0] not in uniqueList: uniqueList.append(row[0]) return uniqueList finally: if row: del row if rows: del rows def getTable(tname): mypath=inWorkspace arcpy.env.workspace =mypath if isshp==True: if arcpy.Exists(mypath+os.sep+tname+".dbf"): return mypath+os.sep+tname+".dbf" else: if arcpy.Exists(mypath+os.sep+tname): return mypath+os.sep+tname return None def getFeatureclass(featurename): mypath=inWorkspace arcpy.env.workspace =mypath if isshp==True: if arcpy.Exists(mypath+os.sep+featurename+".shp"): return mypath+os.sep+featurename+".shp" else: if arcpy.Exists(mypath+os.sep+featurename): return mypath+os.sep+featurename datasets = arcpy.ListDatasets("", "Feature") for dataset in datasets: curpath=mypath+os.sep+dataset if arcpy.Exists(curpath+os.sep+featurename): return curpath+os.sep+featurename return None def updateFieldalias(TableName,FieldName,alias):#更新字段别名 desc = arcpy.Describe(TableName) FieldName=FieldName.upper() for field in desc.fields: if field.Name.upper() ==FieldName: if field.aliasName!=alias: arcpy.AlterField_management(TableName, FieldName, new_field_alias=alias) #修改别名 #field.aliasName=alias arcpy.AddMessage(u"modify'table={2},{0}={1}'".format(FieldName,alias,TableName)) break def FieldExists(TableName,FieldName):#字段是否存在 desc = arcpy.Describe(TableName) for field in desc.fields: if field.Name.upper() ==FieldName.upper(): return True break return False def CreateTable(ftable): num=getCount(ftable) initProgress("create",num) inField=["表名","表英文","类型"] rows = arcpy.da.SearchCursor(ftable,inField) try: for row in rows: step() ptype=row[2] etable=row[1] ctable=row[0] arcpy.AddMessage("etable="+etable+",ctable="+ctable) if ptype==0: table= getTable(etable) if table==None: arcpy.CreateTable_management(inWorkspace, etable) if not isshp:#数据库有别名 arcpy.AlterAliasName(inWorkspace+os.sep+etable, ctable) #修改中文 elif ptype<4: inFeature=getFeatureclass(etable) geometry_type="POINT" if ptype==2: geometry_type="POLYLINE" elif ptype==3: geometry_type="POLYGON" if inFeature==None: #sr = arcpy.Describe("JFB").spatialReference arcpy.CreateFeatureclass_management(inWorkspace, etable, geometry_type, spatial_reference=sr) #arcpy.CreateFeatureclass_management(inWorkspace, etable, geometry_type) #template="#",has_m="DISABLED",has_z="DISABLED",spatial_reference=sr) if not isshp:#数据库有别名 arcpy.AlterAliasName(inWorkspace+os.sep+etable, ctable) #修改中文 finally: freeProgress() if row: del row if rows: del rows def getFieldType(Fieldtype,Fieldlen):#返回标准的字段类型 Fieldtype=Fieldtype.upper() if Fieldtype=="INT" and Fieldlen<5: return "SmallInteger" elif Fieldtype=="INT": return "Integer" elif Fieldtype=="INTEGER": return "Integer" elif Fieldtype=="DOUBLE": return "Double" elif Fieldtype=="FLOAT": return "Double" elif Fieldtype=="STRING": return "String" elif Fieldtype=="CHAR": return "String" else: return "Date" def addoneField(fieldtable,sql,layername,inFeature): rows = arcpy.da.SearchCursor(fieldtable,["字段英文","字段名","字段类型","长度","小数位"],sql,sql_clause=(None,"order by ordid")) try: for row in rows: eField=row[0] cField=row[1] Fieldtype=row[2] Fieldlen=row[3] fieldPrecision=row[4] if not (FieldExists(inFeature,eField)): FType=getFieldType(Fieldtype,Fieldlen) try: if FType.upper()=="Double".upper(): arcpy.AddField_management(inFeature, eField, field_type="DOUBLE",field_precision=Fieldlen,field_scale=fieldPrecision,field_alias=cField) else: arcpy.AddField_management(inFeature, eField, FType, "#","#", Fieldlen, cField) except Exception, ErrorDesc: arcpy.AddWarning(u"错误："+str(ErrorDesc)) finally: if row: del row if rows: del rows def addField(layertable,fieldtable): num=getCount(layertable) inField=["表英文"] initProgress("addField",num) rows = arcpy.da.SearchCursor(layertable,inField) try: for row in rows: step() etable=row[0] inFeature=getFeatureclass(etable) if inFeature==None: continue sql="图层名='"+etable+"'" addoneField(fieldtable,sql,etable,inFeature) finally: freeProgress() if row: del row if rows: del rows def Main(): scriptPath = sys.path[0] mdbpath=scriptPath+os.sep+"Convert.mdb" ftable=mdbpath+os.sep+"图层名" fieldtable=mdbpath+os.sep+"字段" CreateTable(ftable) addField(ftable,fieldtable) inWorkspace=arcpy.GetParameterAsText(0) SR=arcpy.GetParameter(1) #坐标系 #isclockwise=arcpy.GetParameter(4) #是否顺时针，软件自动处理，顺时针 sr = arcpy.SpatialReference() sr.loadFromString(SR) XY=arcpy.GetParameter(2) sr.XYTolerance=XY arcpy.env.XYTolerance=str(XY)+ " Meters" outdesc = arcpy.Describe(inWorkspace) isshp=True if outdesc.dataType== "Workspace": isshp=False Main() elif not outdesc.dataType=="Folder":#如FeatureDataset FeatureLayer arcpy.AddError(u"格式错误无法裁剪") else: Main()
14.6 修改面左上角为第一个点
#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 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): #计算夹角 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) ## extent = geometry.extent ## XMin = extent.XMin ## XMax = extent.XMax ## YMin = extent.YMin ## YMax = extent.YMax 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
14.7 锐角检查
#coding=utf8 import arcpy import os from arcpy import env import math def getCount(inFeature):#获得记录数 result = arcpy.GetCount_management(inFeature) count= int(result.getOutput(0)) return count def getXYAngle(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 splitNgeometry(mgeometry): num=mgeometry.count Sumarray = arcpy.Array() parray = arcpy.Array() for i in range(num): pt=mgeometry[i] if pt: parray.add(pt) else:#内边形 Sumarray.add(parray) parray.removeAll() Sumarray.add(parray) return Sumarray def getLineAnglenew(pt1, pt2, pt3):#获得pt1,pt2,pt3夹角 x1=pt1.X y1=pt1.Y x2=pt2.X y2=pt2.Y x3=pt3.X y3=pt3.Y try: l1 = math.sqrt((x2-x3)*(x2-x3)+(y2-y3)*(y2-y3)) l2 = math.sqrt((x1-x3)*(x1-x3)+(y1-y3)*(y1-y3)) l3 = math.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)) cj =(l1*l1+l3*l3-l2*l2)/(2*l1*l3) return math.acos(cj)* 180 / math.pi except: return 180 ########################################### def getLineAngle(pt1, pt2, pt3,mydir): #获得pt1,pt2,pt3夹角，mydir是true顺时针 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) #arcpy.AddMessage("A1："+str(A1)+",A2："+str(A2)) 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 outputPoint(partgeometry,iCur,pfeat,WN,FID): global RJNUM num=partgeometry.count #arcpy.AddMessage("num:"+str(num)) #if shapeType=="Polyline": num=num-1 for k in range(num): #不包括最后一个 if k==0: if shapeType=="Polyline": continue; if k==0: pt1=partgeometry[num-1] else: pt1=partgeometry[k-1] pt2=partgeometry[k] pt3=partgeometry[k+1] a=getLineAnglenew(pt1,pt2,pt3) #arcpy.AddMessage("K="+str(k)+",角度======:"+str(a)) if a<angle and a>-0.1: pfeat = iCur.newRow() pfeat.setValue(outshapeName, pt2) pfeat.setValue(aField,a) pfeat.setValue(PID,FID) iCur.insertRow(pfeat) arcpy.AddMessage("角度:"+str(a)) RJNUM=RJNUM+1 def getlength(mystr): lenTxt = len(mystr) lenTxt_utf8 = len(mystr.decode('utf-8')) size = int((lenTxt_utf8 - lenTxt)/2 + lenTxt) #arcpy.AddMessage("'{0}'最后长度{1},lenTxt={2},lenTxt_utf8={3}".format(mystr,size,lenTxt,lenTxt_utf8)) return size ##### def midFill(sumn,mystr,Fill): n=getlength(mystr) if n>=sumn: return mystr leftn=int((sumn-n)/2) s="" lefts=s.ljust(leftn,Fill) s="" rightn=sumn-n-leftn rights=s.ljust(rightn,Fill) return lefts+mystr+rights def printauthor(toolname): titlestr="" sumn=60 Fill='*' titlestr=titlestr.ljust(sumn,Fill) arcpy.AddMessage(titlestr) arcpy.AddMessage(midFill(sumn,u"欢迎使用:"+toolname,Fill)) mystr=u"本工具闫磊编制,QQ:276529800,电话:18987281928" arcpy.AddMessage(midFill(sumn,mystr,Fill)) mystr=u"使用前请做好数据备份，工具产生的不良后果请自行承担！" arcpy.AddMessage(midFill(sumn,mystr,Fill)) arcpy.AddMessage(titlestr) def main(): count= getCount(inputFeature) if count <= 0: arcpy.AddMessage(inputFeature+u"没有数据") return arcpy.AddField_management(outFeatures,aField ,"DOUBLE") arcpy.AddField_management(outFeatures,PID ,"LONG") arcpy.SetProgressor("step", u"正在检查数据",0,count,1) rows = arcpy.SearchCursor(inputFeature) iCur, pfeat = None, None iCur = arcpy.InsertCursor(outFeatures) i=0 for row in rows: arcpy.SetProgressorPosition() arcpy.SetProgressorLabel(u"正在等待，完成"+str(round(1.0*i*100/count,2))+"%...") geometry = row.getValue(shapeName) FID=row.getValue(OIDField) part_count = geometry.partCount #有几部分 #arcpy.AddMessage("FID:"+str(FID)+",part_count:"+str(part_count)) for j in range(part_count): partgeometry=geometry.getPart(j) #arcpy.AddMessage("=========================") SplitArray=splitNgeometry(partgeometry) #考虑有内部的面 #arcpy.AddMessage("++++") #SpliArray=splitNgeometry(partgeometry) N=SplitArray.count #arcpy.AddMessage("nnnnnnnnnnnnnnnn="+str(N)) for k in range(N): Splitgeometry=SplitArray[k] outputPoint(Splitgeometry,iCur,pfeat,k,FID) i=i+1 arcpy.ResetProgressor() if (RJNUM>0): arcpy.AddMessage(u"有："+str(RJNUM)+"个,小于"+str(angle)+"度锐角") else: arcpy.AddMessage(u"没有小于"+str(angle)+"度锐角") if iCur: del iCur if pfeat: del pfeat env.overwriteOutput = True inputFeature = arcpy.GetParameterAsText(0) angle = arcpy.GetParameter(1) #角度 outFeatures = arcpy.GetParameterAsText(2) printauthor(u"锐角检查") desc = arcpy.Describe(inputFeature) shapeType=desc.shapeType shapeName = desc.shapeFieldName OIDField=desc.OIDFieldName sr = desc.spatialReference outPath, outFC = os.path.split(outFeatures) arcpy.CreateFeatureclass_management(outPath, outFC, "POINT", "","","",sr) outdesc = arcpy.Describe(outFeatures) outshapeName = outdesc.shapeFieldName aField="angle" PID="PID" RJNUM=0 #锐角格式 main()
14 .8 两面层按面积的叠加面积最大赋属性
#coding=utf8 import arcpy import sys def initProgress(hint,num): arcpy.SetProgressor("step", u"正在"+hint,0,num,1) def step(): arcpy.SetProgressorLabel(u"正在进行....") arcpy.SetProgressorPosition() def freeProgress(): arcpy.ResetProgressor() ########## def getCount(inFeature): result = arcpy.GetCount_management(inFeature) count= int(result.getOutput(0)) return count ##### def delsame(jfb_Select): desc = arcpy.Describe(jfb_Select) shapeName = desc.ShapeFieldName arcpy.DeleteIdentical_management(jfb_Select, [shapeName]) def getOIDField(jfb_Select): desc = arcpy.Describe(jfb_Select) OIDField=desc.OIDFieldName return OIDField ####### def AddLayer(mxd,inFeature): df=arcpy.mapping.ListDataFrames(mxd)[0] addLayer = arcpy.mapping.Layer(inFeature) arcpy.mapping.AddLayer(df, addLayer,"TOP") #AUTO_ARRANGE，"BOTTOM"，TOP def FieldExists(TableName,FieldName): desc = arcpy.Describe(TableName) FieldName=FieldName.upper() for field in desc.fields: if field.Name.upper() ==FieldName: return True break return False def getField(TableName,FieldName): desc = arcpy.Describe(TableName) FieldName=FieldName.upper() for field in desc.fields: if field.Name.upper() ==FieldName: return field break return None def getrowbyFID(FieldNames,FID): with arcpy.da.SearchCursor(byFeature, FieldNames,oidFieldNameby+"="+str(FID)) as cursor: for row in cursor: return row def Main(): count=getCount(inFeature) if count <= 0: arcpy.AddMessage(u""+inFeature+"没有数据") return arcpy.Select_analysis(inFeature,outFeature) #ORIG_FID YL_FID1="YL_FID1" if not FieldExists(inFeature,YL_FID1): arcpy.AddField_management(inFeature,YL_FID1 ,"LONG") FidFieldName=getOIDField(outFeature) arcpy.CalculateField_management(inFeature,YL_FID1, '!'+FidFieldName+'!', "PYTHON") FieldNames= FieldList.split(";") for FieldName in FieldNames: if not FieldExists(outFeature,FieldName): inField=getField(byFeature,FieldName) if inField: arcpy.AddField_management(outFeature,FieldName ,inField.type,inField.precision ,inField.scale,inField.length,inField.AliasName) num=len(FieldNames) for i in range(num-1,-1,-1):#0需要循环的 FieldName=FieldNames[i] #arcpy.AddMessage(u"FieldName"+FieldName+",i="+str(i)) pField=getField(outFeature,FieldName) if not pField.editable or pField.type=="OID" or pField.type=="Geometry": #arcpy.AddMessage(u"FieldName"+FieldName+",只读,删除") FieldNames.pop(i) #移除这个 FieldNames.append("OID@") YL_FID2="YL_FID2" if not FieldExists(byFeature,YL_FID2): arcpy.AddField_management(byFeature,YL_FID2 ,"LONG") FidFieldName=getOIDField(byFeature) arcpy.CalculateField_management(byFeature,YL_FID2, '!'+FidFieldName+'!', "PYTHON") mytemp="in_memory/YL999888" mysort="in_memory/YL999sort" arcpy.AddMessage("TabulateIntersection================") arcpy.TabulateIntersection_analysis(outFeature,YL_FID1,byFeature,mytemp,YL_FID2) arcpy.AddMessage("TabulateIntersection") arcpy.Sort_management(mytemp,mysort,YL_FID1+" ASCENDING;PERCENTAGE DESCENDING") arcpy.AddMessage("Sort_management") arcpy.DeleteIdentical_management(mysort, [YL_FID1]) arcpy.AddMessage("DeleteIdentical") initProgress("update",num) updatecursor=arcpy.da.UpdateCursor(outFeature, FieldNames) scursor=arcpy.da.SearchCursor(mysort, [YL_FID1,YL_FID2,"PERCENTAGE"]) try: srow=scursor.next() Fieldnum=len(FieldNames) k=1 for row in updatecursor: step() FID1=srow[0] FID2=srow[1] Scale=srow[2] #arcpy.AddMessage("========{0},{1},{2}".format(FID1,FID2,Scale)) if Scale>minScale: row2= getrowbyFID(FieldNames,FID2) for i in range(Fieldnum-1): row[i]=row2[i] updatecursor.updateRow(row) else: arcpy.AddMessage(u""+inFeature+"中FID="+str(FID1)+"最大比例"+str(Scale)+",找不到叠加比例大于"+str(minScale)) if k<=num: srow=scursor.next() k=k+1 if FieldExists(byFeature,YL_FID2): arcpy.DeleteField_management(byFeature,YL_FID2) if FieldExists(outFeature,YL_FID1): arcpy.DeleteField_management(outFeature,YL_FID1) if updatecursor: del updatecursor if scursor: del scursor finally: freeProgress() arcpy.env.overwriteOutput = True inFeature = arcpy.GetParameterAsText(0) #输入要素 byFeature = arcpy.GetParameterAsText(1) #依据的要素 FieldList = arcpy.GetParameterAsText(2) #字段列表 outFeature = arcpy.GetParameterAsText(3) #输出要素 minScale=arcpy.GetParameter(4) ##最小比例，如果是负值，取最大的 oidFieldNameby=getOIDField(byFeature) try: Main() except Exception as e: arcpy.AddError(e.message)