# *.* ###############
# -- ##########################
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)
arcpy 狭长图斑检查融合处理工具
于 2024-04-30 09:14:13 首次发布