C#+AE 栅格表面分析

分手后的思念是犯贱 2023-03-13 02:06 11阅读 0赞

### 文章目录 ###

*  一、栅格表面分析
 *   *  1.1坡度
     *  1.2坡向
     *  1.3曲率
 *  二、方法实现
 *   *  2.1坡度计算
     *  2.2坡向计算
     *  2.3等值线计算
     *  2.4填挖方
     *  2.5山体阴影
     *  2.6曲率计算
     *  2.7可见性分析
 *  三、小结

# 一、栅格表面分析 #

> 栅格表面分析，虽然从名字上来看是针对栅格表面来进行分析，但是在AE中这里的栅格表面一般是指“数字地形曲面”(DEM或者DTM)，所以在AE中也可以把栅格表面分析理解为“数字地形分析”。按照这种思维再来理解AE中的很多概念就容易很多了，对于“数字地形分析”而言，其实质上只有两个目的：地形属性计算和地形特征提取，而我们经常会遇到的地形因子本就很多，例如：坡度、坡向、曲率等等，除此之外就是针对地形要素的特征进行提取，这类的地形要素也有很多，例如山脊线、山谷线、流域、可视性等等。

总的来说呢，了解这些地形因子和要素概念有助于我们接下来学习AE的栅格表面分析这部分内容，而且我一直认为在把书上的程序实现的同时，再去真正的理解了程序之外的这些知识点会使我们真正的掌握程序内的东西\*~\*，下面就是我认为有必要了解的几个地形因子。

## 1.1坡度 ##

> 地形表面任一点的坡度是指过该点的切平面与水平面的夹角。其表现了地形表面在该点的倾斜程度，在数值上等于过该点的地表微分单元的法向量n与z轴的夹角。如下图所示：

![20200509191745633.png][]

从上图可以看出坡度值越低，说明该点(像元)所处的地势越平坦，反之，则表明地势较为陡峭。

## 1.2坡向 ##

> 坡向的定义则是地形表面上一点的切平面的法向量在水平面的投影与该点的正北方向的夹角。在上图中，因为是地理坐标系的关系，所以x轴的正半轴所指的方向是正北方向，除此之外特别要说一下的是，坡向表征的是该点高程值改变量的最大变化方向。

我记得我有一段时间一直以为坡向是该点切平面法向量所指的方向,之后仔细看了定义之后才知道自己搞错了啊\*^\*。其实现在想想也能理解一二，像我们现实生活中经常使用东南西北、1点钟方向等等，这些都还停留在二维平面上，而坡向仍然是沿用了这种二维表述，大概是感觉这样描述更能让人理解事物的大致方位。

## 1.3曲率 ##

> 曲率是对地形表面上的一点扭曲变化程度的定量化度量因子，地面曲率在垂直和水平两个方向上的分量则被称为平面曲率和剖面曲率。其中剖面曲率是对地面坡度沿最大坡降方向地面高程变化率的度量，而平面曲率描述的是地表曲面沿水平方向的弯曲、变化情况，也就是该点所在的地面等高线(等值线)的弯曲程度。

![20200509215028397.png][]

因为自己实在是没有绘图的天赋，所以就截取了书上的图来说明了一下这部分内容\*~\*。

# 二、方法实现 #

原始数据图像：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70]

## 2.1坡度计算 ##

代码如下：  
RasterSurfaceAnalysis.cs:

public partial class RasterSurfaceAnalysis : DevExpress.XtraEditors.XtraForm
        { 
            public RasterSurfaceAnalysis()
            { 
                InitializeComponent();
            }
    
            public AxMapControl axMapControl {  get; set; }
            IMap pMap;
    
            private ISurfaceOp2 surfaceOp;       //表面分析操作对象
            private IGeoDataset inputDataset;   //输入数据集
            private IGeoDataset outputDataset;   //输出数据集
            private object Missing = Type.Missing;
    
            private esriGeoAnalysisSlopeEnum slopeEnum;     //输出单位的形式
    
            //窗体生成事件
            private void RasterSurfaceAnalysis_Load(object sender, EventArgs e)
            { 
                surfaceOp = new RasterSurfaceOpClass();
                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 pLayer = GetLayerByName(pMap, comboBox_InputDataset.Text);
                IRasterLayer rasLayer = pLayer as IRasterLayer;
                IRaster raster = rasLayer.Raster;
                inputDataset = raster as IGeoDataset;
            }
    
            /// <summary>
            /// 坡度计算
            /// </summary>
            /// <param name="sender"></param>
            /// <param name="e"></param>
            private void simpleButton_Slope_Click(object sender, EventArgs e)
            { 
                outputDataset = surfaceOp.Slope(inputDataset, slopeEnum, ref Missing);
                SaveAndShowRasterDataset(textEdit_Output.Text);
            }
            //设置输出单位的形式
            private void comboBoxEdit_SlopeUnit_SelectedIndexChanged(object sender, EventArgs e)
            { 
                int index = comboBoxEdit_SlopeUnit.SelectedIndex;
                switch (index)
                { 
                    case 0:
                        slopeEnum = esriGeoAnalysisSlopeEnum.esriGeoAnalysisSlopeDegrees;
                        break;
                    case 1:
                        slopeEnum = esriGeoAnalysisSlopeEnum.esriGeoAnalysisSlopePercentrise;
                        break;
                }
            }
    
            //保存并显示栅格数据集
            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 simpleButton_Output_Click(object sender, EventArgs e)
            { 
                SaveFileDialog flg = new SaveFileDialog();
                flg.Title = "保存路径";
                flg.Filter = "TIFF(.tif)|*.tif|ShpFile(.shp)|*.shp";
                flg.ShowDialog();
    
                textEdit_Output.Text = flg.FileName;
            }
        }

MainForm.cs:

//栅格表面分析
           private void barButtonItem_RasterSurfacceAnalysis_ItemClick(object sender, DevExpress.XtraBars.ItemClickEventArgs e)
           { 
               RasterSurfaceAnalysis surfaceAnalysis = new RasterSurfaceAnalysis();
               surfaceAnalysis.axMapControl = axMapControl_MapView;
               surfaceAnalysis.ShowDialog();
           }

注：输出单位的形式如下：  
![在这里插入图片描述][20200510085953705.png]

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 1]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 2]

## 2.2坡向计算 ##

代码如下：

//坡向计算
            private void simpleButton_Aspect_Click(object sender, EventArgs e)
            { 
                outputDataset= surfaceOp.Aspect(inputDataset);
                SaveAndShowRasterDataset(textEdit_Output.Text);
            }

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 3]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 4]

## 2.3等值线计算 ##

> 在AE中，等值线应理解为“等高线”更为合适，也就是指将栅格表面上相邻且等高程的点相连接而形成的线。

代码如下：

//等值线计算
            private void simpleButton_Contour_Click(object sender, EventArgs e)
            { 
                double interval = Convert.ToDouble(textEdit_Interval.Text.Trim());
                outputDataset = surfaceOp.Contour(inputDataset, interval,ref Missing);
                SaveAndShowFeatureDataset(textEdit_Output.Text);
            }
    
            //保存并显示矢量数据集
            private void SaveAndShowFeatureDataset(string filePath)
            { 
                string directoryPath = System.IO.Path.GetDirectoryName(filePath);
                string[] fileName = System.IO.Path.GetFileName(filePath).Split('.');
    
                IWorkspaceFactory wf = new ShapefileWorkspaceFactoryClass();		//要注意一下这里，因为是要保存矢量文件，所以应该创建矢量工作空间工厂对象
                IWorkspace ws = wf.OpenFromFile(directoryPath, 0) as IWorkspace;
                IConversionOp converop = new RasterConversionOpClass();
                converop.ToFeatureData(outputDataset, ESRI.ArcGIS.Geometry.esriGeometryType.esriGeometryPolyline,
                    ws, fileName[0]);
    
                IFeatureClass featureClass = outputDataset as IFeatureClass;
                IFeatureLayer featLayer = new FeatureLayerClass();
                featLayer.FeatureClass = featureClass;
                featLayer.Name = fileName[0];          //设置图层名字
                axMapControl.AddLayer(featLayer);
                axMapControl.Refresh();
            }

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 5]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 6]

## 2.4填挖方 ##

> 填挖操作是一个通过添加或移除表面物质来修改地表高程的过程。通过给定两个不同的表面，填挖函数会生成一个栅格数据，显示地面物质增加、减少或者不变的区域。负的体积值表明该区域已被填充，正的体积值表明该区域已发生移除。

实现代码如下：

//填挖计算变量
            private IGeoDataset inputAfterDataset;
             /// <summary>
            /// 填挖计算
            /// </summary>
            /// <param name="sender"></param>
            /// <param name="e"></param>
            private void simpleButton_CutFill_Click(object sender, EventArgs e)
            { 
                ILayer pLayer = GetLayerByName(pMap, comboBoxEdit_InputAfterDataset.Text);
                IRasterLayer rasLayer = pLayer as IRasterLayer;
                IRaster raster = rasLayer.Raster;
                inputAfterDataset = raster as IGeoDataset;
    
                outputDataset=surfaceOp.CutFill(inputDataset, inputAfterDataset, ref Missing);
                SaveAndShowRasterDataset(textEdit_Output.Text);
            }

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 7]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 8]  
注：在这里使用IConvensionOp接口保存栅格数据时，发生了异常，说是保存的格式有问题，在将原来的.tif格式更换为.img格式之后，程序正常运行。

## 2.5山体阴影 ##

> 山体阴影操作是指通过考虑照明源的角度和阴影，然后再根据表面栅格来晕染地形地貌的过程。

代码如下：

/// <summary>
            /// 山体阴影
            /// </summary>
            /// <param name="sender"></param>
            /// <param name="e"></param>
            private void simpleButton_HillShade_Click(object sender, EventArgs e)
            { 
                double azimuth = Convert.ToDouble(textEdit_Azimuth.Text.Trim());
                double altitude = Convert.ToDouble(textEdit_Altitude.Text.Trim());
                outputDataset = surfaceOp.HillShade(inputDataset, azimuth, altitude, true, ref Missing);
                SaveAndShowRasterDataset(textEdit_Output.Text);
            }

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 9]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 10]

## 2.6曲率计算 ##

> 曲率是表面的二阶导数，或者可称为坡度的坡度。

实现代码：

/// <summary>
            /// 曲率计算
            /// </summary>
            /// <param name="sender"></param>
            /// <param name="e"></param>
            private void simpleButton_Curvature_Click(object sender, EventArgs e)
            { 
                outputDataset = surfaceOp.Curvature(inputDataset, true, true, ref Missing);
                SaveAndShowRasterDataset(textEdit_Output.Text);
            }

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 11]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 12]

## 2.7可见性分析 ##

实现代码：

//可视性分析
            private IGeoDataset observeDataset;
            private esriGeoAnalysisVisibilityEnum visibilityEnum;
    		/// <summary>
            /// 可视性分析
            /// </summary>
            /// <param name="sender"></param>
            /// <param name="e"></param>
            private void simpleButton_Visibility_Click(object sender, EventArgs e)
            { 
                outputDataset = surfaceOp.Visibility(inputDataset, observeDataset, visibilityEnum);
                SaveAndShowRasterDataset(textEdit_Output.Text);
            }
    
            //设置观测点数据
            private void comboBoxEdit_ObservePoint_SelectedIndexChanged(object sender, EventArgs e)
            { 
                ILayer layer = GetLayerByName(pMap,comboBoxEdit_ObservePoint.Text);
                IFeatureLayer featLayer = layer as IFeatureLayer;
                observeDataset = featLayer.FeatureClass as IGeoDataset;
            }
    
            //分析类型
            private void comboBoxEdit_VisibilityType_SelectedIndexChanged(object sender, EventArgs e)
            { 
                int index = comboBoxEdit_SlopeUnit.SelectedIndex;
                switch (index)
                { 
                    case 0:
                        visibilityEnum = esriGeoAnalysisVisibilityEnum.esriGeoAnalysisVisibilityFrequency;
                        break;
                    case 1:
                        visibilityEnum = esriGeoAnalysisVisibilityEnum.esriGeoAnalysisVisibilityObservers;
                        break;
                }
            }

实现效果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 13]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 14]  
![在这里插入图片描述][202005102025524.png]

# 三、小结 #

整个过程虽然相对简单，但是我也学到了下面两点：

> 1、怎样使用IConversionOp接口去保存栅格数据或者是矢量数据。  
> 2、在AE中的Raster对象和FeatureClass对象才是真正的栅格数据或矢量数据的内存对象，这两个对象会直接控制对栅格数据(或矢量数据)的读取，而他们相应的Dataset对象则是对栅格数据(或矢量数据)的抽象，主要是指磁盘上的数据。

> 参考资料：《地理信息系统-汤国安》《ArcgisEngine地理信息系统开发教程》

[20200509191745633.png]: /images/20230312/6a82affb2ed7418b8e596ccd177b6577.png
[20200509215028397.png]: /images/20230312/e2e30a7396d745fab3759799ac23f8ec.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70]: /images/20230312/7e79b3d0dd0b409093d9f6f5aabe5b5e.png
[20200510085953705.png]: /images/20230312/5ace8cc13a494c5089a912e247704521.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 1]: /images/20230312/45c5576630714fd2b633f52184a11cd2.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 2]: /images/20230312/c9cefbb5ffcc480389536a401fd01fc2.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 3]: /images/20230312/3249418e534843d284c5d6e2cb71eeda.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 4]: /images/20230312/374c2ec9c3e4444cbb9815da344119cd.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 5]: /images/20230312/ac52088ce2bc40f48fae63adbe96d6a2.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 6]: /images/20230312/031225e4da6c4e0897cf33bd48339e47.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 7]: /images/20230312/d65330a6c6254967b242f973296650fa.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 8]: /images/20230312/c6e3d909f8134031a87e08603621e3f3.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 9]: /images/20230312/d62ecaa45eab4c629a23ffaead5a0145.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 10]: /images/20230312/52fe6886c39042c0b94515ac572f1091.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 11]: /images/20230312/e856c7e263fc42b09eda418c7dd22000.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 12]: /images/20230312/c6a0499e534340f8a284d8c2f38a85f2.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 13]: /images/20230312/d329e21d3f6e4b3783c688794dadd7ad.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2RheXVoYWl0YW5nMQ_size_16_color_FFFFFF_t_70 14]: /images/20230312/08e0a86974b44fabb6c503a5b862e521.png
[202005102025524.png]: /images/20230312/660395c313b44775865b46f488a30e83.png