ArcGIS消除叠加操作生成的碎图斑或狭长面

     在ArcGIS中叠加操作会带来碎图斑或狭长面，而这些碎图斑或狭长面又会影响后续的判断。所以就需要消除。

     例如：

具体操作有以下方式：

1、最简单的的方法就是使用ArcGIS提供的“消除工具”。在"数据管理工具"——“制图综合”——“消除”。

2、使用Python脚本

# -*-coding:gbk-*-

import arcpy
from arcpy import env
import time

# 融合指定条件的图斑

try:
    source_gdb_path = "F:/GIS测试数据/测试.gdb"  # 原始图层工作空间
    in_features = "POLYGON"   # 输入图层
    expression = "SHAPE_AREA <50 or SHAPE_Area/ SHAPE_Length<0.2"  # 融合限制条件

    env.workspace = source_gdb_path
    time_begin = time.time()
    print "消除指定条件下的要素{0}".format(expression)

    tempLayer = "block_layer"
    # SearchCursor
    fcount = len([feature[0] for feature in arcpy.da.SearchCursor(in_features, "SHAPE_AREA", expression)])
    while (fcount > 0):
        print('fcount:{0}  Time : {1} s'.format(fcount, time.time() - time_begin))

        tempLayer = "{0}_b".format(tempLayer)
        # MakeFeatureLayer
        arcpy.MakeFeatureLayer_management(in_features, tempLayer)

        in_features = "{0}_E".format(in_features)
        if arcpy.Exists(in_features):
            arcpy.Delete_management(in_features)
        # SelectLayerByAttribute to define feature to be eliminated
        arcpy.SelectLayerByAttribute_management(tempLayer, "NEW_SELECTION", expression)

        # Eliminate
        arcpy.Eliminate_management(tempLayer, in_features, "AREA")
        fcount2 = len([feature[0] for feature in arcpy.da.SearchCursor(in_features, "SHAPE_AREA", expression)])
        if fcount2 == fcount:
            break
        else:
            fcount = fcount2

    print "Total Time {0}".format(time.time() - time_begin)
except Exception, e:
    print e.message.decode("utf_8")

3、AE+C#实现

/// <summary>
/// 消除碎图斑
/// </summary>
/// <param name="pFeatureClass">要素类</param>
/// <param name="whereClause">消除条件："SHAPE_AREA <50 or SHAPE_AREA/SHAPE_Length<0.2"</param>
public static void Eliminate(IFeatureClass pFeatureClass, string whereClause)
{
    Dictionary<int, IGeometry> dicOidsGeo = new Dictionary<int, IGeometry>();
    ISpatialFilter pSpatialFilter = new SpatialFilterClass { WhereClause = whereClause, SubFields = $"SHAPE,{pFeatureClass.OIDFieldName}" };
    IFeatureCursor pFeatureCursor = pFeatureClass.Search(pSpatialFilter, true);
    IFeature pFeature = null;
    while ((pFeature = pFeatureCursor.NextFeature()) != null)
    {
        (pFeature.ShapeCopy as ITopologicalOperator).Simplify();
        dicOidsGeo[pFeature.OID] = pFeature.ShapeCopy;
    }
    System.Runtime.InteropServices.Marshal.ReleaseComObject(pFeatureCursor);
    int count = dicOidsGeo.Count;
    while (count > 0)
    {
        Console.WriteLine($"个数：{count}，时间：{DateTime.Now}");
        int count2 = ExcudeEliminate(dicOidsGeo, pFeatureClass);
        if (count2 == count)
            break;
        else
            count = count2;
    }
}

/// <summary>
/// 消除碎图斑，融到相邻面积最大面
/// </summary>
/// <param name="dicOidsGeo">几何集合</param>
/// <param name="pFeatureClass">源要素类</param>
/// <returns>未消除的个数</returns>
public static int ExcudeEliminate(Dictionary<int, IGeometry> dicOidsGeo, IFeatureClass pFeatureClass)
{
    Func<IFeature, IGeometry, IGeometry> func = (pFeature1, ge) =>
    {
        ITopologicalOperator pTopolo = pFeature1.ShapeCopy as ITopologicalOperator;
        IGeometry geo = pTopolo.Union(ge);
        (geo as ITopologicalOperator).Simplify();
        return geo;
    };

    List<int> lisOids = new List<int>();  //记录已经融合的要素ID

    IFeatureCursor pFeatureCursor = null;
    IFeature pFeature = null;
    ISpatialFilter pSpatialFilter = new SpatialFilterClass { SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects, GeometryField = pFeatureClass.ShapeFieldName, PostfixClause = "Order by SHAPE_AREA desc" };
    foreach (var ge in dicOidsGeo)
    {
        pSpatialFilter.Geometry = ge.Value;
        pFeatureCursor = pFeatureClass.Update(pSpatialFilter, false);
        pFeature = pFeatureCursor.NextFeature();
        if (pFeature != null)
        {
            pFeature.Shape = func(pFeature, ge.Value);
            pFeatureCursor.UpdateFeature(pFeature);

            lisOids.Add(ge.Key);
        }
    }
    System.Runtime.InteropServices.Marshal.ReleaseComObject(pFeatureCursor);

    if (lisOids.Count == 0)
        return dicOidsGeo.Count;
    var str = string.Empty;
    //删除
    lisOids.ForEach(x =>
    {
        str += $",{x}";
        dicOidsGeo.Remove(x);
    });
    str = str.Substring(1);
    IQueryFilter pQueryFilter = new QueryFilterClass();
    pQueryFilter.WhereClause = $"{pFeatureClass.OIDFieldName} in ({str})";  //str超过1000时需要分组
    pFeatureCursor = pFeatureClass.Update(pQueryFilter, false);
    while((pFeature= pFeatureCursor.NextFeature())!= null)
    {
        pFeatureCursor.DeleteFeature();
    }
    System.Runtime.InteropServices.Marshal.ReleaseComObject(pFeatureCursor);
    return dicOidsGeo.Count;
}