C#+AE 空间插值

继续开始ArcgisEngine的学习。

文章目录

• 一、空间插值
• 二、插值方法代码实现
• • 2.1反距离权重法
  • 2.2克里金法
  • 2.3样条函数法
  • 2.4趋势面法
  • 2.5自然邻域法
• 三、小结

一、空间插值

空间数据插值是我们进行数据外推的基本方法，其过程是利用数学曲面来模拟实际地形表面，从而得到我们想要的外推数据。从理论上讲，任何复杂的曲面都可以用多项式进行逼近，但是一个地区常常包含着各种复杂的地形曲面，简单的曲面并不能很好的表达这些地形曲面。所以我们一般在对空间数据进行插值时，都会采用“分而治之”的思想来更好的拟合数学曲面，进而能够更好的表达地形曲面。

在ArcgisEngine中，主要实现了下面这几个空间插值的方法：

1、反距离权重法：这种方法其实从字面上也可以很容易理解，就是将插值点与样本点的距离作为权重来进行插值估计，只不过是这个距离权重是“距离越短，权重越大”。

2、克里金法：该方法是根据相邻变量的值，利用半变异函数(半方差)来揭示区域变量的内在联系进而来预测空间变量的值。

3、样条函数法：该方法是利用最小化表面总曲率的属性函数来估计值，从而生成恰好经过输入点的平滑表面，也就是样条曲面。而所谓的样条曲面可以假想为将一张具有弹性的薄板固定在各个采样点上，而其他的地方自由弯曲的情况。

4、趋势面法：该方法是针对大量离散的点信息，从整体插值出发，进行趋势渐变特征分析的较为简单的方法。其一般是采用多项式来进行回归分析，通过全局多项式插值法来拟合数学曲面。但是因为该方法是整体内插，所以并不常用。

5、自然邻域法：该方法是先找到距离查询点最近的输入样本子集，并根据区域的大小对这些样本再运用权重来进行插值。

二、插值方法代码实现

以下插值方法的实验原始数据图像为：

2.1反距离权重法

新建一个窗体：

SpatialInterpolation.cs:

public partial class SpatialInterpolation : DevExpress.XtraEditors.XtraForm
    {
        public SpatialInterpolation()
        {
            InitializeComponent();
        }

        public AxMapControl axMapControl { get; set; }
        private IMap pMap = null;

        private IRasterAnalysisEnvironment envi;    //创建栅格分析环境对象
        private IInterpolationOp2 interpolation;    //创建空间插值对象
        private IGeoDataset inputDataset;   //输入数据集
        private IGeoDataset outputDataset;   //输出数据集
        private double cellSize = 500;  //像元的大小
        private object cellSizeObj;
        private object Missing = Type.Missing;
        private IRasterRadius rasterRadius; //搜索半径
        IFeatureClass featureClass;
        private object processExtent;   //处理范围
        private IFeatureClassDescriptor featClassDes;   //创建一个要素类描述器对象，其主要是用于控制要素类对象

        //窗体加载事件
        private void SpatialInterpolation_Load(object sender, EventArgs e)
        {
            envi = new RasterInterpolationOpClass();    //实例化栅格插值对象
            pMap = axMapControl.Map;
            if (pMap == null)
                return;
            for (int i = 0; i < pMap.LayerCount; i++)
            {
                comboBox_InputDataset.Properties.Items.Add(pMap.Layer[i].Name);
            }
        }

        //设置输入要素集
        private void comboBox_InputDataset_SelectedIndexChanged(object sender, EventArgs e)
        {
            ILayer currentLayer = GetLayerByName(pMap, comboBox_InputDataset.Text);
            IFeatureLayer featureLayer = currentLayer as IFeatureLayer;
            featureClass = featureLayer.FeatureClass;
            for (int i = 0; i < featureClass.Fields.FieldCount; i++)
            {
                comboBoxEdit_ZValue.Properties.Items.Add(featureClass.Fields.Field[i].Name);
            }
            processExtent = currentLayer;
            envi.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, processExtent, Missing);    //设置空间处理的范围
        }

        //设置Z值字段
        private void comboBoxEdit_ZValue_SelectedIndexChanged(object sender, EventArgs e)
        {
            try
            {
                featClassDes = new FeatureClassDescriptorClass();
                featClassDes.Create(featureClass, null, comboBoxEdit_ZValue.Text);  //创建一个要素类描述器
                inputDataset = featClassDes as IGeoDataset;
            }
            catch
            {
                MessageBox.Show("输入的图层有误！");
            }

        }
        private void simpleButton_IDW_Click(object sender, EventArgs e)
        {
            if (pMap == null)
                return;
            string inputData = comboBox_InputDataset.Text;
            double power = Convert.ToDouble(textEdit_Power.Text.Trim());
            //设置输出栅格大小
            cellSize = Convert.ToDouble(textEdit_CellSize.Text.Trim());
            cellSizeObj = cellSize;
            envi.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, cellSizeObj); 

            //执行IDW方法
            interpolation = envi as IInterpolationOp2;
            outputDataset = interpolation.IDW(inputDataset, power, rasterRadius, Missing);

            //保存输出数据集
            string filePath = textEdit_Output.Text;
            SaveAndShowRasterDataset(filePath);
            MessageBox.Show("执行成功！");
            this.Close();
        }

       //保存并显示栅格数据集
        private void SaveAndShowRasterDataset(string filePath)
        {
            string directoryPath = System.IO.Path.GetDirectoryName(filePath);
            string fileName = System.IO.Path.GetFileName(filePath);

            IWorkspaceFactory wf = new RasterWorkspaceFactoryClass();
            IWorkspace ws = wf.OpenFromFile(directoryPath, 0) as IWorkspace;
            IConversionOp converop = new RasterConversionOpClass();
            converop.ToRasterDataset(outputDataset, "TIFF", ws, fileName);

            IRasterLayer rlayer = new RasterLayerClass();
            IRaster raster = new Raster();
            raster = outputDataset as IRaster;
            rlayer.CreateFromRaster(raster);       //使用raster对象创建一个rasterLayer对象
            rlayer.Name = fileName;    //设置图层名字
            axMapControl.AddLayer(rlayer);
            axMapControl.Refresh();
        }

        //通过图层名获取图层
        private ILayer GetLayerByName(IMap pMap, string layerName)
        {
            ILayer pLayer = null;
            ILayer tempLayer = null;
            for (int i = 0; i < pMap.LayerCount; i++)
            {
                tempLayer = pMap.Layer[i];
                if (tempLayer.Name.ToUpper() == layerName.ToUpper())
                {
                    pLayer = tempLayer;
                    break;
                }
            }
            return pLayer;
        }
        //搜索半径
        private void comboBoxEdit_SearchRadius_SelectedIndexChanged(object sender, EventArgs e)
        {
            rasterRadius = new RasterRadiusClass(); //实例化半径对象
            int index = comboBoxEdit_SearchRadius.SelectedIndex;
            switch (index)
            {
                case 0:
                    rasterRadius.SetFixed(2500, Missing);
                    break;
                case 1:
                    rasterRadius.SetVariable(12, Missing);
                    break;
            }
        }

        /*private void ShowResult(IGeoDataset geoDataset, string layerName)
        {
            IRasterLayer rlayer = new RasterLayerClass();
            IRaster raster = new Raster();
            raster = geoDataset as IRaster;
            rlayer.CreateFromRaster(raster);       //使用raster对象创建一个rasterLayer对象
            rlayer.Name = layerName;    //设置图层名字
            axMapControl.AddLayer(rlayer);
            axMapControl.Refresh();
        }*/

        //输出路径设置
        private void simpleButton_Output_Click(object sender, EventArgs e)
        {
            SaveFileDialog flg = new SaveFileDialog();
            flg.Title = "保存文件";
            flg.Filter = "TIFF(.tif)|*.tif";
            flg.ShowDialog();

            textEdit_Output.Text = flg.FileName;
        }
    }
           

MainForm:

//空间插值
     private void barButtonItem_SpatialInterpolation_ItemClick(object sender, DevExpress.XtraBars.ItemClickEventArgs e)
     {
         SpatialInterpolation interOp = new SpatialInterpolation();
         interOp.axMapControl = axMapControl_MapView;
         interOp.ShowDialog();
     }
           

实现效果：

注：我使用了IConversionOp接口来进行栅格数据的保存，保存的格式为tif文件，如果想要保存成其他文件可以参考下面的详细描述。

2.2克里金法

窗体：

SpatialInterpolation.cs:

//Krige变量
        private esriGeoAnalysisSemiVariogramEnum semiEnum;
 	 	//Krige法
        private void simpleButton_Krige_Click(object sender, EventArgs e)
        {
            //设置输出栅格大小
            cellSize = Convert.ToDouble(textEdit_CellSize.Text.Trim());
            cellSizeObj = cellSize;
            envi.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, cellSizeObj); 

            interpolation = envi as IInterpolationOp2;
            outputDataset= interpolation.Krige(inputDataset, semiEnum, rasterRadius, true, ref Missing);

            SaveAndShowRasterDataset(textEdit_Output.Text);
            MessageBox.Show("执行成功！");
            this.Close();
        }

        //半方差类型
        private void comboBoxEdit_SemiVariogram_SelectedIndexChanged(object sender, EventArgs e)
        {
            int index = comboBoxEdit_SemiVariogram.SelectedIndex;
            switch (index)
            {
                case 0:
                    semiEnum = esriGeoAnalysisSemiVariogramEnum.esriGeoAnalysisNoneVariogram;
                    break;
                case 1:
                    semiEnum = esriGeoAnalysisSemiVariogramEnum.esriGeoAnalysisSphericalSemiVariogram;
                    break;
                case 2:
                    semiEnum = esriGeoAnalysisSemiVariogramEnum.esriGeoAnalysisCircularSemiVariogram;
                    break;
                case 3:
                    semiEnum = esriGeoAnalysisSemiVariogramEnum.esriGeoAnalysisExponentialSemiVariogram;
                    break;
                case 4:
                    semiEnum = esriGeoAnalysisSemiVariogramEnum.esriGeoAnalysisGaussianSemiVariogram;
                    break;
                case 5:
                    semiEnum = esriGeoAnalysisSemiVariogramEnum.esriGeoAnalysisLinearSemiVariogram;
                    break;
            }
        }
           

Krige的半方差参数类型：

实现效果：

2.3样条函数法

代码如下：

//样条函数变量
        private esriGeoAnalysisSplineEnum splintEnum;
        //样条函数
        private void simpleButton_Spline_Click(object sender, EventArgs e)
        {
            //设置输出栅格大小
            cellSize = Convert.ToDouble(textEdit_CellSize.Text.Trim());
            cellSizeObj = cellSize;
            envi.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, cellSizeObj);

            interpolation = envi as IInterpolationOp2;
            outputDataset= interpolation.Spline(inputDataset, splintEnum, ref Missing, ref Missing);

            SaveAndShowRasterDataset(textEdit_Output.Text);
            MessageBox.Show("执行成功！");
            this.Close();
        }
        private void comboBoxEdit_Spline_SelectedIndexChanged(object sender, EventArgs e)
        {
            int index = comboBoxEdit_Spline.SelectedIndex;
            switch (index)
            {
                case 0:
                    splintEnum = esriGeoAnalysisSplineEnum.esriGeoAnalysisRegularizedSpline;
                    break;
                case 1:
                    splintEnum = esriGeoAnalysisSplineEnum.esriGeoAnalysisTensionSpline;
                    break;
            }
        }
           

实现效果：

2.4趋势面法

代码如下：

//Trend变量
        private esriGeoAnalysisTrendEnum trendEnum;
	    /// <summary>
        /// 趋势面法
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void simpleButton_Trend_Click(object sender, EventArgs e)
        {
            int power = Convert.ToInt32(textEdit_Power.Text.Trim());    //多项式的阶数
            //设置输出栅格大小
            cellSize = Convert.ToDouble(textEdit_CellSize.Text.Trim());
            cellSizeObj = cellSize;
            envi.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, cellSizeObj);

            interpolation = envi as IInterpolationOp2;
            outputDataset = interpolation.Trend(inputDataset, trendEnum, power);

            SaveAndShowRasterDataset(textEdit_Output.Text);
            MessageBox.Show("执行成功！");
            this.Close();
        }
        //设置回归类型
        private void comboBoxEdit_TrendType_SelectedIndexChanged(object sender, EventArgs e)
        {
            int index = comboBoxEdit_TrendType.SelectedIndex;
            switch (index)
            {
                case 0:
                    trendEnum = esriGeoAnalysisTrendEnum.esriGeoAnalysisLinearTrend;
                    break;
                case 1:
                    trendEnum = esriGeoAnalysisTrendEnum.esriGeoAnalysisLogisticTrend;
                    break;
            }
        }
           

实现效果：

2.5自然邻域法

代码如下：

private void simpleButton_NatureNeighbor_Click(object sender, EventArgs e)
        {
            //设置输出栅格大小
            cellSize = Convert.ToDouble(textEdit_CellSize.Text.Trim());
            cellSizeObj = cellSize;
            envi.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, cellSizeObj);

            interpolation = envi as IInterpolationOp2;
            outputDataset = interpolation.NaturalNeighbor(inputDataset);

            SaveAndShowRasterDataset(textEdit_Output.Text);
            MessageBox.Show("执行成功！");
            this.Close();
        }
           

实现效果：

三、小结

整个流程走下来，其实我们的工作就是在想办法来创建每个方法所需要的参数对象，然后再借助IInterpolationOp2接口来执行相应的方法即可。但是在这个过程中，我自我感觉是应该尽可能的去理解这些参数更深层的意义，这样不仅有助于我们可以更好的把握整个程序，而且还有助于我们去更好的理解这部分的一些专业知识*~*。

参考资料：《ArcgisEngine地理信息系统开发教程》