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    }

上海市的处理结果如下：