需求
有一面状地类图斑要素类,需要检测出所有位于要素面范围内的孔洞,根据其面积大小进行不同的处理(比如,大于50亩需要人工复核,小于50亩报错)。
思路
- 第一种思路:逐个遍历,检测面中是否包含内洞。具体方式是将其转为GeometryCollection对象,判断是否具有面积为负值的多边形,再判断面积。此种方式执行效率较低,速度慢,尤其是图斑数量多而孔洞少的时候。
- 第二种思路:利用GP工具,先面为线,再线转面,再与原图层求几何差集,结果即为所有孔洞。遍历孔洞,与原图层进行空间查询,取得包含该洞的多边形后,根据孔洞面积判断处理方式。对图斑数量多而孔洞数据较少时,此方法适用。
实现
先利用Python脚本处理数据,得到孔洞要素类,并预先计算好面积字段值(单位:亩):
# -*-coding:gbk-*-
# 提取成果库城镇村等用地图层,取得“空洞”,转为面要素,并计算
# 以“亩”为单位的面积。
#
# 备注:调用时,遍历TEMP_HOLE图层取得空洞要素,判断面积,与
# 成果库城镇村等用地图层进行空间相交分析,取得空洞所在要素信息
import arcpy
import os
import sys
import time
result_gdb_path = sys.argv[1] # 成果库GDB路径
result_name = sys.argv[2] # 成果库中待分析图层名(城镇村等用地)
temp_gdb_path = sys.argv[3] # 临时GDB路径
# result_gdb_path = r"E:/Code/ibuild_gdsdzj/build/CacheData/2001H2019440607.gdb"
# result_name = r"CZCDYD"
# temp_gdb_path = r"E:/Temp/SD_Temp.gdb"
if __name__ == '__main__':
try:
time_begin = time.time()
temp_result_name = "{0}_R".format(result_name) # 临时库中成果库要素类名
temp_polyline_name = "{0}_L".format(result_name) # 临时库中成果库要素类转多段线要素类名
temp_polygon_name = "{0}_G".format(result_name) # 临时库中成果库要素类转多边形要素类名
temp_hole_name = "TEMP_HOLE" # 临时库中处理结果(空洞要素面,带单位为“亩”的面积)要素类名
arcpy.env.workspace = temp_gdb_path
# 删除已存在的图层
for name in [temp_result_name, temp_polyline_name, temp_polygon_name, temp_hole_name]:
if arcpy.Exists(name):
print "Delete {0}".format(name)
arcpy.Delete_management(name)
# 拷贝数据到临时库
print "Copy {0} to {1}".format(result_name, temp_result_name)
arcpy.Copy_management(os.path.join(result_gdb_path, result_name), os.path.join(temp_gdb_path, temp_result_name))
# 要素转线
print "Feature {0} to Line {1}".format(temp_result_name, temp_polyline_name)
arcpy.FeatureToLine_management(temp_result_name, temp_polyline_name)
# 要素转面
print "Feature {0} to Polygon {1}".format(temp_polyline_name, temp_polygon_name)
arcpy.FeatureToPolygon_management(temp_polyline_name, temp_polygon_name)
# 两个面取几何差集
print "SymDiff_analysis {0},{1} to {2}".format(temp_result_name, temp_polygon_name, temp_hole_name)
arcpy.SymDiff_analysis(temp_result_name, temp_polygon_name, temp_hole_name, "ONLY_FID")
# 添加“亩”字段并计算
field_hole_area = "AREA_MU"
print "Add Field {0} of {1}".format(field_hole_area, temp_hole_name)
arcpy.AddField_management(temp_hole_name, field_hole_area, "DOUBLE")
print "Calculate Field {0} of {1}".format(field_hole_area, temp_hole_name)
arcpy.CalculateField_management(temp_hole_name, field_hole_area, "!Shape_area!*0.0015", "PYTHON_9.3")
# 删除中间图层(只保留结果图层)
for name in [temp_result_name, temp_polyline_name, temp_polygon_name]:
if arcpy.Exists(name):
print "Delete {0}".format(name)
arcpy.Delete_management(name)
print "Total Time {0}".format(time.time() - time_begin)
print "1"
except Exception, e:
print e.message
最后, 在C#写的ArcGIS Engine代码中,遍历结果图层,取得孔洞进行相交查询,即可获得逐个有孔洞图斑的信息,并预以处理。
//用孔洞逐个去相交查询城镇村等用地图层要素
//孔洞面积大于阈值(50亩)时需人工复核表,否则报错
using (var comReleaser = new ComReleaser())
{
var pFeatureClassContainer = pSouthDataset as IFeatureClass;
var pDataset = pFeatureClassContainer as IDataset;
//生成孔洞
var pFeatureClassHole = GenerateHoleFeatureClass((pFeatureClassContainer as IDataset).Workspace.PathName, pDataset.Name);
comReleaser.ManageLifetime(pFeatureClassHole);
//查询孔洞,取面积(亩数)
var sFieldHoleArea = "AREA_MU";
var pSearchCursorHole = pFeatureClassHole.Search(new QueryFilterClass() { SubFields = $"{sFieldHoleArea},Shape" }, true);
comReleaser.ManageLifetime(pSearchCursorHole);
var nIndexHoleArea = pFeatureClassHole.FindField(sFieldHoleArea);
IFeature pFeatureHole = null;
comReleaser.ManageLifetime(pFeatureHole);
//以孔洞为几何图形,相交查询城镇村等用地范围层,取标识码和OID
var sFieldContainerGid = "BSM";
var nIndexContainerGid = pFeatureClassContainer.FindField(sFieldContainerGid);
IFeature pFeatureContainer = null;
comReleaser.ManageLifetime(pFeatureContainer);
IFeatureCursor pFeatureCursorContainer = null;
comReleaser.ManageLifetime(pFeatureCursorContainer);
var pSpatialFilterContainer = new SpatialFilterClass()
{
SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects
};
comReleaser.ManageLifetime(pSpatialFilterContainer);
var setBigContainerOid = new HashSet<int>(); //已经输出过人工复核信息的OID集合
var setSmallContainerOid = new HashSet<int>(); //已经输出过错误信息的OID集合
IGeometryCollection pGeometryCollection = null;
//遍历孔洞,查询城镇村等用地,输出人工复核或错误信息
while (null != (pFeatureHole = pSearchCursorHole.NextFeature()))
{
var dAreaHole = Math.Abs((pFeatureHole.Shape as IArea).Area);
var bBigHole = (double)pFeatureHole.Value[nIndexHoleArea] > dArea; //大于指定阈值,需要人工复核,否则认定为错误
pSpatialFilterContainer.Geometry = pFeatureHole.Shape;
pFeatureCursorContainer = pFeatureClassContainer.Search(pSpatialFilterContainer, true);
while (null != (pFeatureContainer = pFeatureCursorContainer.NextFeature()))
{
pGeometryCollection = pFeatureContainer.Shape as IGeometryCollection;
var nGeoCount = pGeometryCollection.GeometryCount;
if (2 > nGeoCount)
continue;
var bContainer = false;
for (int i = 0; i < nGeoCount; i++)
{
if (!(pGeometryCollection.Geometry[i] is IArea pArea))
continue;
var dAreaGeometry = Math.Abs(pArea.Area);
if (Math.Abs(dAreaGeometry - dAreaHole) > 1e-4)
continue;
bContainer = true;
break;
}
if (!bContainer)
continue;
var nOid = pFeatureContainer.OID;
if (bBigHole) //大于阈值时,已报过人工复核消息的跳过,否则添加到集合
{
if (setBigContainerOid.Contains(nOid))
continue;
else
setBigContainerOid.Add(nOid);
}
else //小于阈值时,已报过错误消息的跳过,否则添加到集合
{
if (setSmallContainerOid.Contains(nOid))
continue;
else
setSmallContainerOid.Add(nOid);
}
var sGid = pFeatureContainer.Value[nIndexContainerGid] + string.Empty;
if (bBigHole)
{
//大于50亩的处理
}
else
{
//小于50亩的处理
}
}
}
}