ArcGIS栅格图像重分类C#代码实现

最近在做热点商圈分析，在分析出来的结果要对其进行重分类，通过不断尝试，终于实现了以不同方式（NaturalBreaks、EqualInterval、GeometricalInterval、Quantile）进行重分类的功能，现在将其整理分享：

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.DataSourcesGDB;
using ESRI.ArcGIS.esriSystem;
using System.IO;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.Carto;
using ESRI.ArcGIS.GeoAnalyst;

namespace Intensity.Reclassify
{
    public partial class frmReclassify : Form
    {
        //define global varibale
        IWorkspaceFactory pWorkspaceFactory = new FileGDBWorkspaceFactoryClass();
        IWorkspace pWorkspace = null;

        IEnumDataset pEnumDataset = null;
        IEnumDatasetName pEnumDatasetName = null;

        IClassifyGEN mClassify;
        int mClassesNo;//declare reclassify count
        object varNewValues, varNewCounts;//declare variable to store pixel values


        string strFileGDBPath = "";
        string strReclassifyResultPath = "";

        public frmReclassify()
        {
            InitializeComponent();
        }
        /// <summary>
        /// form load
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void frmReclassify_Load(object sender, EventArgs e)
        {
            cmbClassifyMethod.Items.Add("NaturalBreaks");
            cmbClassifyMethod.Items.Add( "EqualInterval");
            cmbClassifyMethod.Items.Add("GeometricalInterval");
            cmbClassifyMethod.Items.Add("Quantile");
        }
        /// <summary>
        /// Select Raster Dataset GDB file
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void btnSelectRasterGDB_Click(object sender, EventArgs e)
        {
            FolderBrowserDialog dlg = new FolderBrowserDialog();
            dlg.Description = "打开GDB文件";
            if (DialogResult.OK == dlg.ShowDialog())
            {
                if (Directory.Exists(dlg.SelectedPath))
                {
                    strFileGDBPath = dlg.SelectedPath;
                    txtFileGDBPath.Text = strFileGDBPath;
                }
                else
                {
                    return;
                }
                pWorkspace = pWorkspaceFactory.OpenFromFile(strFileGDBPath, 0);
                pEnumDataset = pWorkspace.get_Datasets(esriDatasetType.esriDTRasterDataset);

                //pRasterWorkspace = pWorkspace as IRasterWorkspace;

                if (pEnumDataset == null)
                {
                    MessageBox.Show("文件地理库中不包含栅格数据，请检查！");
                    return;
                }
            }
        }
        /// <summary>
        /// Select Classify Method
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void cmbClassifyMethod_SelectedIndexChanged(object sender, EventArgs e)
        {
            string strMethod = this.cmbClassifyMethod.Text.Trim();
            if (strMethod == "") return;
            switch (strMethod)
            {
                case "EqualInterval":
                    {
                        mClassify = new EqualIntervalClass();
                        break;
                    }
                case "GeometricalInterval":
                    {
                        mClassify = new GeometricalIntervalClass();
                        break;
                    }
                case "NaturalBreaks":
                    {
                        mClassify = new NaturalBreaksClass();
                        break;
                    }
                case "Quantile":
                    {
                        mClassify = new QuantileClass();
                        break;
                    }
                default:
                    {
                        mClassify = new EqualIntervalClass();
                        break;
                    }
            }
        }
        /// <summary>
        /// select Reclassify GDB path
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void btnSelectResultPath_Click(object sender, EventArgs e)
        {
            FolderBrowserDialog dlg = new FolderBrowserDialog();
            dlg.Description = "打开GDB文件";
            if (DialogResult.OK == dlg.ShowDialog())
            {
                if (Directory.Exists(dlg.SelectedPath))
                {
                    strReclassifyResultPath = dlg.SelectedPath;
                    this.txtReclassifyPath.Text = strReclassifyResultPath;
                }
                else
                {
                    return;
                }
                
            }
        }
        /// <summary>
        /// Start to Reclassify the dataset
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void btnOK_Click(object sender, EventArgs e)
        {
            if (this.txtFileGDBPath.Text.Trim() == "" || this.cmbClassifyMethod.Text.Trim() == "" || txtClassifyCount.Text.Trim() == "")
            {
                MessageBox.Show("请先设置重分类条件！");
                return;
            }

            int iReclassifyCount = Convert.ToInt16(this.txtClassifyCount.Text.Trim());
            string strReclassifyMethod = this.cmbClassifyMethod.Text.Trim();

            pEnumDataset.Reset();
            IRasterDataset pRasterDataset = new RasterDatasetClass();
            pRasterDataset = pEnumDataset.Next() as IRasterDataset;
            frmProgress frmPro = new frmProgress();
            frmPro.Show();
            frmPro.TopMost = true;
            frmPro.SetLabel1("正在计算,请等待......");
            frmPro.SetProgressRefresh();
            frmPro.AutoSize = true;
            frmPro.Refresh();

            Application.DoEvents();
            while (pRasterDataset != null)
            {
                IRasterLayer pRasterLayer = new RasterLayerClass();
                pRasterLayer.CreateFromDataset(pRasterDataset);
               

                frmPro.SetLabel1("正在生成" + pRasterLayer.Name + "的Reclassify图,请等待......");
                try
                {
                    ReclassifyRasterLayerAndSave(pRasterLayer, iReclassifyCount, strReclassifyMethod, strReclassifyResultPath, "Reclassify_" + pRasterLayer.Name);

                }
                catch
                {
                    MessageBox.Show(pRasterLayer.Name + "生成Reclassify不成功，请检查!");
                    if (frmPro != null)
                    {
                        frmPro.Close();
                        frmPro.Dispose();
                    }
                }
                pRasterDataset = pEnumDataset.Next() as IRasterDataset;


            }
            frmPro.SetLabel1("Reclassify Complete!");
        }
        /// <summary>
        /// exit the program
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void btnExit_Click(object sender, EventArgs e)
        {
            this.Close();
        }

        

    #region Customize method
            private void ReclassifyRasterLayerAndSave(IRasterLayer pRasterLayer, int pClassesNo, string pMethodName, string strOutputPath, string strRasterName)
            {
                try
                {
                    mClassesNo = Convert.ToInt16(this.txtClassifyCount.Text.Trim());

                    IRaster pRaster = pRasterLayer.Raster;
                    IRasterProps rasterProps = (IRasterProps)pRaster;
                    //定义字典集合存放像素值和统计频数

                    Dictionary<double, long> ValueFrequence = new Dictionary<double, long>();

                    //设置栅格数据起始点  
                    IPnt pBlockSize = new Pnt();
                    pBlockSize.SetCoords(rasterProps.Width, rasterProps.Height);
                    //选取整个范围  
                    IPixelBlock pPixelBlock = pRaster.CreatePixelBlock(pBlockSize);
                    //左上点坐标  
                    IPnt tlp = new Pnt();
                    tlp.SetCoords(0, 0);
                    //读入栅格  
                    IRasterBandCollection pRasterBands = pRaster as IRasterBandCollection;
                    IRasterBand pRasterBand = pRasterBands.Item(0);
                    IRawPixels pRawRixels = pRasterBands.Item(0) as IRawPixels;
                    pRawRixels.Read(tlp, pPixelBlock);


                    //将PixBlock的值组成数组  
                    System.Array pSafeArray = pPixelBlock.get_SafeArray(0) as System.Array;

                    for (int y = 0; y < rasterProps.Height; y++)
                    {
                        for (int x = 0; x < rasterProps.Width; x++)
                        {
                            double value = Convert.ToDouble(Convert.ToDouble(pSafeArray.GetValue(x, y)).ToString("##.0000000000"));


                            if (!ValueFrequence.ContainsKey(value))
                            {
                                ValueFrequence.Add(value, 1);
                            }
                            else if (ValueFrequence.ContainsKey(value))
                            {
                                ValueFrequence[value] = ValueFrequence[value] + 1;
                            }


                        }
                    }

                    DataTable dt = initialDatatable();

                    foreach (double dValue in ValueFrequence.Keys)
                    {
                        DataRow pRow = dt.NewRow();
                        pRow["PixelValue"] = dValue;
                        pRow["PixelFrequence"] = long.Parse(ValueFrequence[dValue].ToString());
                        dt.Rows.Add(pRow);
                    }
		    //排序，很重要
                    dt.DefaultView.Sort = "PixelValue asc,PixelFrequence";

                    DataTable ddt = dt.DefaultView.ToTable();
                    int iCount = ddt.Rows.Count;
                    double[] dValues = new double[iCount];
                    long[] lFrequency = new long[iCount];
                    for (int ii = 0; ii < iCount; ii++)
                    {
                        dValues[ii] = Convert.ToDouble(ddt.Rows[ii]["PixelValue"].ToString());
                        lFrequency[ii] = long.Parse(ddt.Rows[ii]["PixelFrequence"].ToString());
                    }

                    varNewValues = dValues as object;
                    varNewCounts = lFrequency as object;


                    if (mClassify == null) return;

                    //mClassify 为IClassifyGEN接口对象，
                    //通过其子类进行初始化实现mClassify = new NaturalBreaksClass();
                    //进行分类（varNewValues--像素值数组；varNewCounts--对应值的个数数组；
                    //pClassesNo--分类级数）
                    mClassify.Classify(varNewValues, varNewCounts, ref pClassesNo);

                    //分类完成,获得断点
                    double[] ClassNum = (double[])mClassify.ClassBreaks;

                    IRasterDescriptor pRD = new RasterDescriptorClass();
                    IReclassOp pReclassOp = new RasterReclassOpClass();
                    IGeoDataset pGeodataset = pRasterLayer.Raster as IGeoDataset;

                    //IRasterAnalysisEnvironment pEnv = pReclassOp as IRasterAnalysisEnvironment;
                    //pEnv.OutWorkspace = pWorkspace;
                    //object objSnap = null;
                    //object objExtent = pGeodataset.Extent;
                    //pEnv.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, ref objExtent, ref objSnap);
                    //pEnv.OutSpatialReference = pGeodataset.SpatialReference;

                    IRasterLayer pRLayer = new RasterLayerClass();

                    INumberRemap pNumRemap = new NumberRemapClass();
                    for (int i = 0; i < mClassesNo; i++)
                    {
                        pNumRemap.MapRange(ClassNum[i], ClassNum[i + 1], i + 1);
                    }
                    IRemap pRemap = pNumRemap as IRemap;

                    IRaster pOutRaster = pReclassOp.ReclassByRemap(pGeodataset, pRemap, false) as IRaster;

                    pRLayer.CreateFromRaster(pOutRaster);

                    IRasterBandCollection pRasterBandCol = pOutRaster as IRasterBandCollection;
                    IRasterBand pBand = pRasterBandCol.Item(0);
                    IRasterDataset pOutRasterDataset = pBand.RasterDataset;

                    IWorkspaceFactory pp = new FileGDBWorkspaceFactoryClass();
                    IWorkspace ppWorkspace = pp.OpenFromFile(strOutputPath, 0);

                    ISaveAs pSaveAs = pOutRasterDataset as ISaveAs;
                    pSaveAs.SaveAs(strRasterName, ppWorkspace, "GDB");
                }
                catch (Exception ex)
                {
                    MessageBox.Show(ex.Message);
                    return;
                }

            }

            public IRasterWorkspaceEx OpenFileGDBWorkspace(string sPath)
            {
                IWorkspaceFactory workspaceFactory = new FileGDBWorkspaceFactoryClass();
                IWorkspace workspace = workspaceFactory.OpenFromFile(sPath, 0);
                IRasterWorkspaceEx rasterWorkspaceEx = (IRasterWorkspaceEx)workspace;
                return rasterWorkspaceEx;
            }

            private DataTable initialDatatable()
            {
                try
                {
                    DataTable dt = new DataTable();
                    DataColumn col = new DataColumn();
                    col.ColumnName = "PixelValue";
                    col.DataType = System.Type.GetType("System.Double");
                    dt.Columns.Add(col);

                    col = new DataColumn();
                    col.ColumnName = "PixelFrequence";
                    col.DataType = System.Type.GetType("System.Int64");
                    dt.Columns.Add(col);

                    return dt;
                }
                catch
                {
                    return null;
                }
            }
    #endregion
    }

上海市的处理结果如下：