zoukankan      html  css  js  c++  java
  • 14章例子

    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)
  • 相关阅读:
    Flesch Reading Ease(模拟)
    实验一:词法分析设计
    java—容器学习笔记
    [转载]马士兵Java视频教程 —— 学习顺序
    Java的安装过程
    编程之美初赛第一场
    RCC 2014 Warmup (Div. 2)
    ural 1017. Staircases(dp)
    ural 1012. K-based Numbers. Version 2(大数dp)
    ural 1009. K-based Numbers(简单dp)
  • 原文地址:https://www.cnblogs.com/gisoracle/p/13524136.html
Copyright © 2011-2022 走看看