arcpy 狭长图斑检查融合处理工具

# *.* ###############
# -- ##########################
import random
import re
import string
import sys
import arcpy
import os

from arcpy import Point
import json

reload(sys)
sys.setdefaultencoding('utf-8')

intable = arcpy.GetParameterAsText(0)  # 输入图层
outtable = arcpy.GetParameterAsText(1)  # 输出结果
area_limit = arcpy.GetParameterAsText(2)
SD_limit = arcpy.GetParameterAsText(3)
# intable = r'C:\temp\test.gdb\Export_Output1'  # 输入图层
# outtable = r'C:\temp\test.gdb\Export_Output5'  # 输出结果
area_limit = 1
SD_limit = 0.1


def tableconversion_REMOVE(intable, prefix, outtable, IsTrue=True):
    fieldmappings = arcpy.FieldMappings()
    infields = arcpy.ListFields(intable)
    outname = outtable[outtable.rfind('/') + 1:]
    outpath = outtable[:outtable.rfind('/')]
    for infield in infields:
        if infield.required:
            continue
        fieldmap = arcpy.FieldMap()
        fieldmap.addInputField(intable, infield.name)
        outfieldname = re.sub(prefix, '', infield.name, flags=re.IGNORECASE)
        outfield = fieldmap.outputField
        outfield.isNullable = IsTrue
        outfield.name = outfieldname
        fieldmap.outputField = outfield
        fieldmappings.addFieldMap(fieldmap)
        del fieldmap, outfield
    # if arcpy.Exists(outname):
    #     arcpy.Delete_management(outname)
    arcpy.FeatureClassToFeatureClass_conversion(intable, out_path=outpath, out_name=outname,
                                                field_mapping=fieldmappings)


def tableconversion_ADD(intable, prefix, outtable, IsTrue=True):
    fieldmappings = arcpy.FieldMappings()
    infields = arcpy.ListFields(intable)
    outname = outtable[outtable.rfind('/') + 1:]
    outpath = outtable[:outtable.rfind('/')]
    for infield in infields:
        if infield.required:
            continue
        fieldmap = arcpy.FieldMap()
        fieldmap.addInputField(intable, infield.name)
        if prefix.upper() not in infield.name.upper():
            outfieldname = prefix + infield.name
        else:
            outfieldname = infield.name
        outfield = fieldmap.outputField
        # outfield.isNullable = IsTrue
        outfield.name = outfieldname
        fieldmap.outputField = outfield
        fieldmappings.addFieldMap(fieldmap)
        del fieldmap, outfield
    # if arcpy.Exists(outname):
    #     arcpy.Delete_management(outname)
    arcpy.FeatureClassToFeatureClass_conversion(intable, out_path=outpath, out_name=outname,
                                                field_mapping=fieldmappings)


def isFieldExist(fc, fieldname):
    fieldList = arcpy.ListFields(fc)
    for field in fieldList:
        if (field.name.upper() == fieldname.upper()):
            return True
    return False


def AlterFieldsName(inTable, prefix):
    fields = arcpy.ListFields(inTable)
    for field in fields:
        if not field.required:
            arcpy.AlterField_management(inTable, field.name, field.name.replace(prefix, ''))


# 检查是否狭长图�?def checkTBLength(perimeter, area, SD_limit):
    result = 2 * pow(3.14 * area, 0.5) / perimeter
    if result < SD_limit:
        return True
    else:
        return False


def CaculateAngle(pointArray, angle):
    try:
        for i in range(0, len(pointArray)):
            print i, pointArray[i]
            ptL = Point()
            ptM = Point()
            ptR = Point()
            if i < len(pointArray) - 2:
                ptL = pointArray[i]
                ptM = pointArray[i + 1]
                ptR = pointArray[i + 2]
            elif i == len(pointArray) - 2:
                ptL = pointArray[i]
                ptM = pointArray[i + 1]
                ptR = pointArray[0]
            elif i == len(pointArray) - 1:
                ptL = pointArray[i]
                ptM = pointArray[0]
                ptR = pointArray[1]
                # continue
            AML = Point()
            AML.X = ptL.X - ptM.X
            AML.Y = ptL.Y - ptM.Y

            AMR = Point()
            AMR.X = ptR.X - ptM.X
            AMR.Y = ptR.Y - ptM.Y

            ab = (AML.X * AMR.X + AML.Y * AMR.Y)
            M = (math.sqrt(AML.X * AML.X + AML.Y * AML.Y) * math.sqrt(AMR.X * AMR.X + AMR.Y * AMR.Y))
            if M == 0:
                continue
            cosA = ab / M

            # 向量乘积小于0，说明是大于90度，不用判断
            if ab < 0:
                continue
            # 因为0到π/2是递减,所以小于cos给定角度值，那么角度比给定角度值大,就返回了
            # elif cosA > math.cos(angle / 180 * math.pi):
            elif cosA > math.cos(math.radians(angle)):
                return True
            elif cosA <= math.cos(angle / 180 * math.pi):
                continue
        # 遍历所有的顶点了，说明满足要求
        return False
    except Exception as ex:
        print ex.message


try:
    intable = intable.replace('\\', '/')
    outtable = outtable.replace('\\', '/')

    m_intable = intable[intable.rfind('/') + 1:]

    temppath = os.getenv('TEMP')
    # temppath = "c:\\temp"
    gdb = arcpy.CreateFileGDB_management(temppath, "gdb_%s" % "".join(
        random.sample(string.uppercase + string.digits, 5)) + '.gdb')
    tmppath = str(gdb) + '\\' + m_intable
    arcpy.Copy_management(intable, tmppath)

    fieldList = arcpy.ListFields(tmppath)
    field_name_List = []
    for field in fieldList:
        field_name_List.append(field.name)
    if 'Shape_Area' not in field_name_List:
        arcpy.AddField_management(tmppath, 'Shape_Area', 'DOUBLE', field_precision=18, field_scale=6)
        arcpy.CalculateField_management(tmppath, 'Shape_Area', '!shape.area!', 'PYTHON_9.3')
    if 'Shape_Length' not in field_name_List:
        arcpy.AddField_management(tmppath, 'Shape_Length', 'DOUBLE', field_precision=18, field_scale=6)
        arcpy.CalculateField_management(tmppath, 'Shape_Length', '!shape.length!', 'PYTHON_9.3')

    # edit = arcpy.da.Editor(str(gdb))
    # edit.startEditing(True, True)
    # edit.startOperation()

    mfattribute_table = arcpy.MakeFeatureLayer_management(tmppath)
    mflocation_table = arcpy.MakeFeatureLayer_management(tmppath)
    dic_area = {}

    with arcpy.da.SearchCursor(tmppath, ['shape_area', 'Shape_Length', 'OBJECTID', 'bsm']) as cursor:
        for row in cursor:
            # if row[0] < area_limit or (checkTBLength(row[1], row[0], SD_limit)):  # 小于最小上图面积或是大于最小上图面积且符合细碎图斑
            if row[0] < area_limit:
                flag = True
                mfattribute_table = arcpy.MakeFeatureLayer_management(tmppath)
                arcpy.SelectLayerByAttribute_management(mfattribute_table, selection_type='NEW_SELECTION',
                                                        where_clause='objectid =' + str(row[2]))  # 处理被裁切成细碎图斑的宗地要�?                mflocation_table = arcpy.MakeFeatureLayer_management(tmppath)
                arcpy.SelectLayerByLocation_management(mflocation_table, 'BOUNDARY_TOUCHES', mfattribute_table,
                                                       selection_type='NEW_SELECTION')  # 处理被裁切成细碎图斑的宗地要�?                dic_area = {}

                with arcpy.da.SearchCursor(mflocation_table, ['shape_area', 'objectid', 'shape@', 'BSM']) as c:
                    for r in c:
                        dic_area[r[1]] = [r[0], r[1], arcpy.FromWKT(r[2].WKT), r[3]]

                target_geo = dic_area[row[2]][2]
                targetobjectid = row[2]

                maxinterlength = 0  # 求最长公共边
                dic_maxlenth = {}

                # 优先匹配最大公共边
                for key, value in dic_area.items():
                    if key != row[2] and value[3] == dic_area[row[2]][3]:
                        tmplength = target_geo.intersect(value[2], 2).length
                        if maxinterlength == 0 or maxinterlength < tmplength:
                            dic_maxlenth = {}
                            maxinterlength = tmplength
                            dic_maxlenth[key] = value[0]
                        elif tmplength == maxinterlength:
                            dic_maxlenth[key] = value[0]

                if len(dic_maxlenth.values()) >= 1:  # 多个公共边长度相等的情况下，取面积最大的
                    sortdic = sorted(dic_maxlenth.items(), key=lambda d: d[1], reverse=True)
                else:
                    flag = False
                if flag:
                    isfinish = False
                    for s in sortdic:
                        issucceed = False

                        with arcpy.da.UpdateCursor(tmppath,
                                                   ['shape_area', 'objectid', 'shape@'],
                                                   where_clause='objectid=' + str(s[0])) as uc:

                            for ur in uc:
                                checkresult = True
                                maxobjectid = s[0]
                                max_geo = dic_area[maxobjectid][2]
                                newunion = max_geo.union(target_geo)
                                json_obj = json.loads(newunion.JSON)
                                for rings in json_obj['rings']:
                                    pointArray = []
                                    for dd in rings:
                                        point = Point()
                                        point.X = dd[0]
                                        point.Y = dd[1]
                                        pointArray.append(point)
                                    if CaculateAngle(pointArray, 10):
                                        checkresult = False
                                        break
                                if checkresult:
                                    ur[2] = newunion
                                    uc.updateRow(ur)
                                    issucceed = True
                                    break
                                else:
                                    continue
                        if issucceed:  # 融合的图形未出现小于10度的夹角，则视为成功，则删除狭长图斑
                            with arcpy.da.UpdateCursor(tmppath,
                                                       ['shape_area', 'objectid', 'shape@'],
                                                       where_clause='objectid=' + str(targetobjectid)) as uc:
                                for ur in uc:
                                    uc.deleteRow()
                            isfinish = True
                            break
                    # 最长公共边无法融合
                    if not isfinish:
                        dic_maxlenth = {}
                        for key, value in dic_area.items():
                            if key != row[2] and value[3] == dic_area[row[2]][3]:
                                tmplength = target_geo.intersect(value[2], 2).length
                                dic_maxlenth[key] = tmplength
                        if len(dic_maxlenth.values()) >= 1:  # 按公共边长降序排�?                            sortdic = sorted(dic_maxlenth.items(), key=lambda d: d[1], reverse=True)
                            del (sortdic[0])
                        for s in sortdic:
                            issucceed = False
                            with arcpy.da.UpdateCursor(tmppath,
                                                       ['shape_area', 'objectid', 'shape@'],
                                                       where_clause='objectid=' + str(s[0])) as uc:

                                for ur in uc:
                                    checkresult = True
                                    maxobjectid = s[0]
                                    max_geo = dic_area[maxobjectid][2]
                                    newunion = max_geo.union(target_geo)
                                    json_obj = json.loads(newunion.JSON)
                                    for rings in json_obj['rings']:
                                        pointArray = []
                                        for dd in rings:
                                            point = Point()
                                            point.X = dd[0]
                                            point.Y = dd[1]
                                            pointArray.append(point)
                                        if CaculateAngle(pointArray, 10):
                                            checkresult = False
                                            break
                                    if checkresult:
                                        ur[2] = newunion
                                        uc.updateRow(ur)
                                        issucceed = True
                                        break
                                    else:
                                        continue
                            if issucceed:  # 融合的图形未出现小于10度的夹角，则视为成功，则删除狭长图斑
                                with arcpy.da.UpdateCursor(tmppath,
                                                           ['shape_area', 'objectid', 'shape@'],
                                                           where_clause='objectid=' + str(targetobjectid)) as uc:
                                    for ur in uc:
                                        uc.deleteRow()
                                break

    arcpy.Copy_management(tmppath, outtable)

except Exception as ex:
    arcpy.AddError(ex.message)