基于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;        }

直接上图吧～