使用ArcPy进行图斑内洞检查

需求

有一面状地类图斑要素类，需要检测出所有位于要素面范围内的孔洞，根据其面积大小进行不同的处理（比如，大于50亩需要人工复核，小于50亩报错）。

思路

1. 第一种思路：逐个遍历，检测面中是否包含内洞。具体方式是将其转为GeometryCollection对象，判断是否具有面积为负值的多边形，再判断面积。此种方式执行效率较低，速度慢，尤其是图斑数量多而孔洞少的时候。
2. 第二种思路：利用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亩的处理
			}
		}
	}
}