基于ArcEngine进行地物分类景观指数计算

前两天写了一个对地物分类进行景观评价的小工具，采用ArcEngine实现，在此备忘一下，仅供学习参考.

核心过程如下：根据评价尺度大小分块获取影响，分别对各个分块进行地物统计，计算景观指数，将该得分保存在分块中心点位置.

        // rst        : 原始地物分类影像
        // patchSize  : 分块大小，正方形边长，为该计算单元包含像元个数
        // shpFile    : 结果保存路径，点状要素类，保存为ShpeFile文件
        private void runModel(Raster rst, int patchSize, string shpFile)
        {
            #region 准备工作
            
            // 获取背景值
            dynamic nodata = ((rst as IRasterBandCollection).Item(0) as IRasterProps).NoDataValue;
            appendMessage("\tNodataValue = " + nodata);

            // 初始化分块辅助工具
            var rstProps = rst as IRasterProps;
            var pixelSize = rstProps.MeanCellSize();
            var rstExt = rstProps.Extent;
            var pbHelper = new PixelBlockHelper(rstProps.Width, rstProps.Height,pixelSize.X, pixelSize.Y,
                        rstExt.XMin, rstExt.XMax, rstExt.YMin, rstExt.YMax);
            long pbCount = pbHelper.GetCount(patchSize);
            int pbCols = pbHelper.GetColumns(patchSize);
            int pbRows = pbHelper.GetRows(patchSize);
            appendMessage("\tPatch Number : " + pbRows + "(Row) X " + pbCols + "(Col) = " + pbCount);
            // 初始化分块游标
            IRaster2 rst2 = rst as IRaster2;
            var patchCursor = rst2.CreateCursorEx(new Pnt() { X = patchSize, Y = patchSize });

            // 创建shapefile文件并开启编辑状态
            var scoreName = "pjz";
            var resultShp = CommonFuns.pointShp(shpFile, rstProps.SpatialReference, scoreName);
            IFeatureClassWrite shpWriter = resultShp as IFeatureClassWrite;
            IWorkspaceEdit shpWpsEdit = (shpWriter as IDataset).Workspace as IWorkspaceEdit;
            shpWpsEdit.StartEditing(true);
            shpWpsEdit.StartEditOperation();
            var scoreFieldIndex = resultShp.FindField(scoreName);

            // 进度指示器
            long doneStep = pbCount / 10;                           // 每隔 10% 播报进度
            long doneCount = 0;                                     // 当前进度

            GeoLandEvaluator glEvaluator = new GeoLandEvaluator();  // 初始化地物评价工具
            IFeature resultFeature;                                 // 评价结果要素
            #endregion

            #region 地物景观指数计算
            ESRI.ArcGIS.esriSystem.ISet featureSet = new ESRI.ArcGIS.esriSystem.Set();
            do
            {
                glEvaluator.Reset();                                // 重置评价工具
                var block = patchCursor.PixelBlock as IPixelBlock3;
                dynamic lands = block.get_PixelData(0);             // 获取分快像元值，System.Array
                for (int x = 0; x < block.Width; x++)               // 扫描像元值填充入评价工具，过滤背景值
                {
                    for (int y = 0; y < block.Height; y++)
                    {
                        dynamic landType = lands.GetValue(x, y);
                        if (nodata != landType) glEvaluator.AddLand((int)landType);
                    }
                }
                var score = glEvaluator.Score();                    // 计算得分
                var tl = patchCursor.TopLeft;                       // 计算中心点位置
                var cxy = pbHelper.GetCenter(tl.X, tl.Y, patchSize);

                // 保存计算结果,播报进度
                resultFeature = resultShp.CreateFeature();
                resultFeature.Shape = new ESRI.ArcGIS.Geometry.Point()
                {
                    X = cxy[0],
                    Y = cxy[1]
                };
                resultFeature.set_Value(scoreFieldIndex, score);
                featureSet.Add(resultFeature); 
                doneCount++;
                if (doneCount % doneStep == 0)
                {
                    shpWriter.WriteFeatures(featureSet);
                    featureSet.RemoveAll();
                    appendMessage("\t...  " + (100.0 * doneCount / pbCount) + " %  ...");
                }
            } while (patchCursor.Next());
            shpWriter.WriteFeatures(featureSet);
            featureSet.RemoveAll();

            #endregion

            // 关闭编辑状态，保存数据
            shpWpsEdit.StopEditOperation();
            shpWpsEdit.StopEditing(true);

            appendMessage("<< Finished! ");
        }

地物评价工具工作过程为：统计分块内各种类型的地物面积大小，计算各类地物的百分比 Pi ，采用如下公式计算景观指数得分——

 score = -（P1＊ln(P1)+P2*ln(P2)+...+Pn*ln(Pn)）

using System;
using System.Collections.Generic;

namespace ...
{
    /// <summary>
    /// 地物类别统计评估工具
    /// </summary>
    public class GeoLandEvaluator
    {
        /// <summary>
        /// 地物分类计数器
        /// </summary>
        private IDictionary<int, long> landCounter = null;

        /// <summary>
        /// 地物总数
        /// </summary>
        private long landSum = 0;

        public GeoLandEvaluator()
        {
            landCounter = new Dictionary<int, long>();
            Reset();
        }

        /// <summary>
        /// 重置地物统计评估工具
        /// </summary>
        public void Reset()
        {
            if (landCounter != null) landCounter.Clear(); 
            landSum = 0;
        }

        /// <summary>
        /// 添加分类
        /// </summary>
        public void AddLand(int landType)
        {
            if (landCounter.ContainsKey(landType))
            {
                landCounter[landType]++;
            }
            else
            {
                landCounter.Add(landType, 1);
            }
            landSum++;
        }

        /// <summary>
        /// 添加地物分类
        /// </summary>
        public void AddLand(int[] landTypes)
        {
            foreach (int ldType in landTypes)
            {
                AddLand(ldType);
            }
        }

        /// <summary>
        /// 计算景观指数评价值
        /// </summary>
        /// <returns></returns>
        public double Score()
        {
            if (landCounter == null || landCounter.Count < 1) return 0;
            double score = 0.0;
            foreach (long ldCount in landCounter.Values)
            {
                double per = 1.0 * ldCount / landSum;
                score += per * Math.Log(per, Math.E);
            }
            return score * (-1.0);
        }
    }
}

栅格像元分块辅助工具

using System;

namespace ...
{
    /// <summary>
    /// 栅格像元分块辅助工具
    /// </summary>
    public class PixelBlockHelper
    {
        public double PixelWidth { get; private set; }
        public double PixelHeight { get; private set; }

        public int RstWidth { get; private set; }
        public int RstHeight { get; private set; }

        public double RstXMin { get; private set; }
        public double RstYMin { get; private set; }
        public double RstXMax { get; private set; }
        public double RstYMax { get; private set; }


        public PixelBlockHelper(int rstWidth, int rstHeight,double pixelWidth, double pixelHeight,
            double rstXmin, double rstXmax, double rstYmin, double rstYmax)
        {
            this.RstWidth = rstWidth;
            this.RstHeight = rstHeight;
            this.PixelWidth = pixelWidth;
            this.PixelHeight = pixelHeight;
            this.RstXMin = rstXmin;
            this.RstXMax = rstXmax;
            this.RstYMin = rstYmin;
            this.RstYMax = rstYmax;
        }

        /// <summary>
        /// 计算 分块 中心坐标 ， 计算原点为影像（可视）左上角
        /// </summary>
        /// <param name="xOffset">左上角横轴偏移量</param>
        /// <param name="yOffset">左上角纵轴偏移量</param>
        /// <param name="blockSize">正方形分块大小，边长</param>
        /// <returns>分块中心点 地图坐标</returns>
        public double[] GetCenter(double xOffset, double yOffset, int blockSize)
        {
            double[] xyOffset = new double[] { xOffset, yOffset };
            return GetCenter(xyOffset, blockSize);
        }


        /// <summary>
        /// 计算 分块 中心坐标 ， 计算原点为影像（可视）左上角
        /// </summary>
        /// <param name="xyOffset">左上角偏移量，[x,y]</param>
        /// <param name="blockSize">正方形分块大小，边长</param> 
        /// <returns>分块中心点 地图坐标</returns>
        public double[] GetCenter(double[] xyOffset, int blockSize)
        {
            xyOffset[0] = PixelWidth * (xyOffset[0] + blockSize / 2.0) + RstXMin;
            xyOffset[1] = RstYMax - PixelHeight * (xyOffset[1] + blockSize / 2.0);
            return xyOffset;
        }

        /// <summary>
        /// 计算分块 总数
        /// </summary>
        /// <returns>可划分的大小为 blockSize 的正方形分块总数</returns>
        public long GetCount(int blockSize)
        {
            return (long)GetColumns(blockSize) * GetRows(blockSize);
        }

        /// <summary>
        /// 计算X方向 分块 列数
        /// </summary>
        /// <returns>以blockSize为间隔可划分的列数</returns>
        public int GetColumns(int blockSize)
        {
            return (int)Math.Ceiling((double)(this.RstWidth) / blockSize);
        }

        /// <summary>
        /// 计算 Y方向 分块 行数
        /// </summary>
        /// <returns>以blockSize为间隔可划分的行数</returns>
        public int GetRows(int blockSize)
        {
            return (int)Math.Ceiling((double)(this.RstHeight) / blockSize);
        }
    }
}

创建Shpefile文件核心代码如下：

/// <summary>
        /// 创建 点 Shapefile
        /// </summary>
        /// <param name="shpFile">完整的保存路径</param>
        /// <param name="spatialRef">参考坐标系</param>
        /// <param name="scoreName">属性值字段名称（double类型字段），默认为 pjz</param>
        /// <returns>包含字段类型为double名称为scoreName的shapefile要素类</returns>
        public static IFeatureClass pointShp(string shpFile, ISpatialReference spatialRef, string scoreName = "pjz")
        {
            FileInfo fi = new FileInfo(shpFile);
            var wpsDir = fi.DirectoryName;   // shapefile保存目录
            var shapeFiledName = "Shape";
            var fidFieldName = "FID";
            if (string.IsNullOrEmpty(scoreName)) scoreName = "pjz";
            //
            var factory = new ShapefileWorkspaceFactory() as IWorkspaceFactory;
            var shpWps = (IFeatureWorkspace)factory.OpenFromFile(wpsDir, 0);

            // prepare fileds for the shapefile
            Fields fieldCollection = new Fields();
            IFieldsEdit fieldCollectionEdit = fieldCollection as IFieldsEdit;

            // 设置FID
            Field field = new Field();
            IFieldEdit fieldEdit = field as IFieldEdit;
            fieldEdit.Name_2 = fidFieldName;
            fieldEdit.Type_2 = esriFieldType.esriFieldTypeOID;
            fieldCollectionEdit.AddField(field);

            // 设置几何类型字段
            field = new Field();
            fieldEdit = field as IFieldEdit;
            fieldEdit.Name_2 = shapeFiledName;
            fieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;
            GeometryDef geoDef = new GeometryDef();
            IGeometryDefEdit geoDefEdit = geoDef as IGeometryDefEdit;
            geoDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPoint;
            geoDefEdit.SpatialReference_2 = spatialRef;
            fieldEdit.GeometryDef_2 = geoDef;
            fieldCollectionEdit.AddField(field);

            // 设置景观指数评价值字段
            field = new Field();
            fieldEdit = field as IFieldEdit;
            fieldEdit.Name_2 = scoreName;
            fieldEdit.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldCollectionEdit.AddField(field);

            //
            var shp = shpWps.CreateFeatureClass(fi.Name, fieldCollection, null, null,
                      esriFeatureType.esriFTSimple,shapeFiledName, string.Empty);
            return shp;
        }

直接上图吧～