查找错误的等值线中的高程点
算法经过了两次判断，第一次，高程点与所在三角形的顶点之间的判断（ok, error, unknown），第二次，针对第一次中的unknown，如果将三角形逐层外推，利用高程点高程、高程点所在三角形的顶点高程，以及外推结果的三角形的顶点高程进行比较
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# ErrorPointFinder.py
# Created on: 2016-07-28 19:10:17.00000
# (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------
# Set the necessary product code
# import arcinfo
import json
import time
# Import arcpy module
import arcpy
# Check out any necessary licenses
arcpy.CheckOutExtension("3D")
arcpy.CheckOutExtension("spatial")
# Set Geoprocessing environments
ProjCRSXian80Proj = "PROJCS['Xian_1980_3_Degree_GK_CM_120E',GEOGCS['GCS_Xian_1980',DATUM['D_Xian_1980',SPHEROID['Xian_1980',6378140.0,298.257]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Gauss_Kruger'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',120.0],PARAMETER['Scale_Factor',1.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]"
arcpy.env.outputCoordinateSystem = ProjCRSXian80Proj
arcpy.env.overwriteOutput = True
# Local variables:
TestPolyline = "E:\\data\\testGP\\testService.gdb\\TestData\\Elevation_L"
TestPoint = "E:\\data\\testGP\\testService.gdb\\TestData\\Elevation_P"
testTIN = "E:\\data\\testGP\\testtin"
trianglePolygonSHP = "E:\\data\\testGP\\output\\trianglePolygonTemp.shp"
trianglePolygonLocation = "E:\\data\\testGP\\output\\out.gdb"
trianglePolygonName = "trianglePolygon"
trianglePolygon = trianglePolygonLocation + "\\" + trianglePolygonName
trianglePolyline = "E:\\data\\testGP\\output\\out.gdb\\trianglePolyline"
spatialJoinResult = "E:\\data\\testGP\\output\\out.gdb\\decisionResult"
# Process: Calculate Field
codeBlockMaxZ='''
def MaxZ(shape):
line = shape.getPart(0)
pnt = line.next()
maxValue = float("-inf")
while pnt:
if maxValue < pnt.Z:
maxValue = pnt.Z
pnt = line.next()
return maxValue
'''
codeBlockMinZ='''
def MinZ(shape):
line = shape.getPart(0)
pnt = line.next()
minValue = float("inf")
while pnt:
if minValue > pnt.Z:
minValue = pnt.Z
pnt = line.next()
return minValue
'''
codeBlockMeanZ='''
def MeanZ(shape):
line = shape.getPart(0)
pnt = line.next()
meanValue = 0
count = 0
while pnt:
count += 1
meanValue += pnt.Z
pnt = line.next()
return meanValue/count
'''
codeBlockFindErrorInfo='''
def FindErrorInfo( ZValue , ZValueMin , ZValueMax ):
returnValue = "unknown"
if ZValue and ZValueMin and ZValueMax:
if ZValueMax - ZValueMin < 0.00000001:
returnValue = "unknown"
elif ZValue > ZValueMax or ZValue < ZValueMin:
returnValue = "error"
else:
returnValue = "ok"
return returnValue
'''
#-------------------------------------------
def getFieldMappings():
fieldMappings = arcpy.FieldMappings()
#
fieldMapZ = arcpy.FieldMap()
fieldMapZ.addInputField(TestPoint, "Z")
outputFieldZ = fieldMapZ.outputField
outputFieldZ.name = 'Z'
outputFieldZ.name = 'Z'
fieldMapZ.outputField = outputFieldZ
fieldMappings.addFieldMap(fieldMapZ)
#
fieldMapTriIndex = arcpy.FieldMap()
fieldMapTriIndex.addInputField(trianglePolygon, "Tri_Index")
outputFieldTriIndex = fieldMapTriIndex.outputField
outputFieldTriIndex.name = 'TriIndex'
outputFieldTriIndex.name = 'TriIndex'
fieldMapTriIndex.outputField = outputFieldTriIndex
fieldMappings.addFieldMap(fieldMapTriIndex)
#
fieldMapZValueMax = arcpy.FieldMap()
fieldMapZValueMax.addInputField(trianglePolygon, "ZValueMax")
outputFieldZValueMax = fieldMapZValueMax.outputField
outputFieldZValueMax.name = 'ZValueMax'
outputFieldZValueMax.name = 'ZValueMax'
fieldMapZValueMax.outputField = outputFieldZValueMax
fieldMappings.addFieldMap(fieldMapZValueMax)
#
fieldMapZValueMin = arcpy.FieldMap()
fieldMapZValueMin.addInputField(trianglePolygon, "ZValueMin")
outputFieldZValueMin = fieldMapZValueMin.outputField
outputFieldZValueMin.name = 'ZValueMin'
outputFieldZValueMin.name = 'ZValueMin'
fieldMapZValueMin.outputField = outputFieldZValueMin
fieldMappings.addFieldMap(fieldMapZValueMin)
#
fieldMapZValueMean = arcpy.FieldMap()
fieldMapZValueMean.addInputField(trianglePolygon, "ZValueMean")
outputFieldZValueMean = fieldMapZValueMean.outputField
outputFieldZValueMean.name = 'ZValueMean'
outputFieldZValueMean.name = 'ZValueMean'
fieldMapZValueMean.outputField = outputFieldZValueMean
fieldMappings.addFieldMap(fieldMapZValueMean)
#
fieldMapFeatureID = arcpy.FieldMap()
fieldMapFeatureID.addInputField(trianglePolygon, "FeatureID")
outputFieldFeatureID = fieldMapFeatureID.outputField
outputFieldFeatureID.name = 'FeatureID'
outputFieldFeatureID.name = 'FeatureID'
fieldMapFeatureID.outputField = outputFieldFeatureID
fieldMappings.addFieldMap(fieldMapFeatureID)
return fieldMappings
def preprocessData():
# Process: Create TIN
sourceShape = TestPolyline + " Z Hard_Line Z"
arcpy.CreateTin_3d(testTIN, ProjCRSXian80Proj, sourceShape, "CONSTRAINED_DELAUNAY")
print("finished CreateTin_3d")
# Process: TIN Triangle
arcpy.TinTriangle_3d(testTIN, trianglePolygonSHP, "DEGREE", "1")
print("finished TinTriangle_3d")
# Process: Feature Class to Feature Class
fieldMappings = arcpy.FieldMappings()
fieldMappings.addTable(trianglePolygonSHP)
arcpy.FeatureClassToFeatureClass_conversion(trianglePolygonSHP, trianglePolygonLocation, trianglePolygonName, "",fieldMappings,"")
print("finished FeatureClassToFeatureClass_conversion")
# Process: Add Field ZValueMax ZValueMin ZValueMean FeatureID
arcpy.AddField_management(trianglePolygon, "ZValueMax", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(trianglePolygon, "ZValueMin", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(trianglePolygon, "ZValueMean", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(trianglePolygon, "FeatureID", "LONG", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
# Process: CalculateField ZValueMax ZValueMin ZValueMean FeatureID
arcpy.CalculateField_management(trianglePolygon, "ZValueMax", "MaxZ(!SHAPE!)", "PYTHON_9.3", codeBlockMaxZ)
arcpy.CalculateField_management(trianglePolygon, "ZValueMin", "MinZ(!shape!)", "PYTHON_9.3", codeBlockMinZ)
arcpy.CalculateField_management(trianglePolygon, "ZValueMean", "MeanZ(!shape!)", "PYTHON_9.3", codeBlockMeanZ)
arcpy.CalculateField_management(trianglePolygon, "FeatureID", "!OBJECTID!", "PYTHON_9.3", "")
print("finished CalculateField_management")
# Process: Polygon To Line
arcpy.PolygonToLine_management(trianglePolygon, trianglePolyline, "IDENTIFY_NEIGHBORS")
print("finished PolygonToLine_management")
# Process: Spatial Join
fieldMappings = getFieldMappings()
arcpy.SpatialJoin_analysis(TestPoint, trianglePolygon, spatialJoinResult, "JOIN_ONE_TO_MANY", "KEEP_ALL", fieldMappings, "INTERSECT", "", "")
print("finished SpatialJoin_analysis")
# Process: Add Field ErrorInfo ErrorInfoRefine
arcpy.AddField_management(spatialJoinResult, "errorInfo", "TEXT", "", "", "50", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(spatialJoinResult, "errorInfoRefine", "TEXT", "", "", "50", "", "NULLABLE", "NON_REQUIRED", "")
# Process: Calculate Field
arcpy.CalculateField_management(spatialJoinResult, "errorInfo", "FindErrorInfo( !Z! , !ZValueMin! , !ZValueMax! )", "PYTHON_9.3", codeBlockFindErrorInfo)
print("finished CalculateField_management")
#----------------------------------------------------------------------------
def getTopologicData():
#
pointInfoList = []
cursor = arcpy.da.SearchCursor(spatialJoinResult, ["TARGET_FID", "Z", "FeatureID", "ZValueMin", "ZValueMax", "ZValueMean", "errorInfo"], "\"Join_Count\" > 0")
for row in cursor:
info = {"TFID": row[0], "Z":row[1], "FID":row[2], "ZVMin":row[3], "ZVMax":row[4], "ZVMean":row[5], "errorInfo":row[6]}
#print(info)
pointInfoList.append(info)
del cursor
polygonInfoList = []
cursor = arcpy.da.SearchCursor(trianglePolygon, ["FeatureID", "ZValueMin", "ZValueMax", "ZValueMean"], "1=1")
for row in cursor:
info = {"FID":row[0], "ZVMin":row[1], "ZVMax":row[2], "ZVMean":row[3]}
#print(info)
polygonInfoList.append(info)
del cursor
polylineInfoList = []
cursor = arcpy.da.SearchCursor(trianglePolyline, ["LEFT_FID", "RIGHT_FID"], "1=1")
for row in cursor:
info = {"Lfid":row[0], "Rfid":row[1]}
#print(info)
polylineInfoList.append(info)
del cursor
'''
f = open("pointInfoList.txt", 'r');
pointInfoList = json.load(f)
#print(pointInfoList)
f.close
#
f = open("polygonInfoList.txt", 'r');
polygonInfoList = json.load(f)
#print(polygonInfoList)
f.close
#
f = open("polylineInfoList.txt", 'r');
polylineInfoList = json.load(f)
#print(polylineInfoList)
f.close
'''
print("finished SearchCursor")
return pointInfoList, polygonInfoList, polylineInfoList
#----------------------------------------------------------------------------
#
def getPolygon(polygonInfoList, FID):
for polygonInfo in polygonInfoList:
if polygonInfo["FID"] == FID:
return polygonInfo
return None
#
def isVisited(visitedPolygon, FID):
flag = False
for polygonInfo in visitedPolygon:
if polygonInfo["FID"] == FID:
flag = True
break
return flag
#
def getNearPolygon(FID, polygonInfoList, polylineInfoList, visitedPolygon):
listResult = []
for polylineInfo in polylineInfoList:
#
if polylineInfo["Lfid"] == FID and polylineInfo["Rfid"] > 0:
flag = False
newFID = polylineInfo["Rfid"]
flag = isVisited(visitedPolygon, newFID)
if flag == False:
polygonInfo = getPolygon(polygonInfoList, newFID)
visitedPolygon.append(polygonInfo)
listResult.append(polygonInfo)
#
if polylineInfo["Rfid"] == FID and polylineInfo["Lfid"] > 0:
flag = False
newFID = polylineInfo["Lfid"]
flag = isVisited(visitedPolygon, newFID)
if flag == False:
polygonInfo = getPolygon(polygonInfoList, newFID)
visitedPolygon.append(polygonInfo)
listResult.append(polygonInfo)
return listResult
def decisionRules(pointValue, polygonInnerValue, outerPolygonInfo):
decisionInfo = "error"
'''
if pointValue < polygonInnerValue and polygonInnerValue < outerPolygonInfo:
decisionInfo = "ok"
elif pointValue > polygonInnerValue and polygonInnerValue > outerPolygonInfo:
decisionInfo = "ok"
'''
# "ZVMin":row[1], "ZVMax":row[2], "ZVMean":row[3]}
if pointValue < polygonInnerValue["ZVMin"] and pointValue > outerPolygonInfo["ZVMin"]:
decisionInfo = "ok"
if pointValue > polygonInnerValue["ZVMax"] and pointValue < outerPolygonInfo["ZVMax"]:
decisionInfo = "ok"
if abs(polygonInnerValue["ZVMin"] - outerPolygonInfo["ZVMin"]) < 0.000001 and abs(polygonInnerValue["ZVMax"] - outerPolygonInfo["ZVMax"]) < 0.000001:
decisionInfo = "unknown"
return decisionInfo
def refineProcessData():
pointInfoList, polygonInfoList, polylineInfoList, = getTopologicData()
refineInfoList = []
for pointInfo in pointInfoList:
if pointInfo["errorInfo"] == "ok" or pointInfo["errorInfo"] == "error":
continue
print(pointInfo)
errorInfo = "unknown"
visitedPolygon=[]
pointValue = pointInfo["Z"]
polygonInnerFID = pointInfo["FID"]
polygonInfoInner = getPolygon(polygonInfoList, polygonInnerFID)
if polygonInfoInner == None:
print("no polygon")
continue
visitedPolygon.append(polygonInfoInner)
#
nearPolygonList = []
nearPolygonList.append(polygonInfoInner)
level = 1
while 1==1:
print("level: " + str(level))
outerPolygonList = []
for currentPolygon in nearPolygonList:
FID = currentPolygon["FID"]
currentNearPolygon = getNearPolygon(FID, polygonInfoList, polylineInfoList, visitedPolygon)
for polygonInfo in currentNearPolygon:
outerPolygonList.append(polygonInfo)
print("level: " + str(level) + " outerPolygonList: " + str(len(outerPolygonList)))
# no nearPolygon
if len(outerPolygonList)<1:
break
#
for outerPolygonInfo in outerPolygonList:
errorInfo = decisionRules(pointValue, polygonInfoInner, outerPolygonInfo)
if errorInfo == "ok":
break
if errorInfo == "error" or errorInfo == "ok":
break
level += 1
print(errorInfo)
refineInfoList.append({"TARGET_FID": pointInfo["TFID"], "errorInfoRefine":errorInfo})
print("finished refineProcessData")
return refineInfoList
#
def updateData(refineInfoList):
cursor = arcpy.UpdateCursor(spatialJoinResult)
for row in cursor:
tfid = row.getValue("TARGET_FID")
for refineInfo in refineInfoList:
if tfid == refineInfo["TARGET_FID"]:
row.setValue("errorInfoRefine", refineInfo["errorInfoRefine"])
cursor.updateRow(row)
del row
del cursor
print("finished UpdateCursor")
#
preprocessData()
refineInfoListRe = refineProcessData()
updateData(refineInfoListRe)
#
print("-------------------finish---------------------")
