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地理資訊系統開發教程》