AE二次开发进行DEM表面积求算

1、.四、算法代码此后为自定义类中的代码1栅格数据提取：采用二维数组对栅格数据进展结果存储，代码为arcengine帮助文档内部搜索出然后修正的public double, ReadWriteRawBlocksIRasterDataset2 rasterDs            /Create a raster.             IRaster2 raster2 = ra

2、sterDs.CreateFullRaster as IRaster2;            /Create a raster cursor with a system-optimized pixel block size by passing a null.            IRasterCursor rasterCursor = raster2.CreateCurs

3、orExnull;            /Use the IRasterEdit interface.            IRasterEdit rasterEdit = raster2 as IRasterEdit;            /Loop throu

4、gh each band and pixel block.            IRasterBandCollection bands = rasterDs as IRasterBandCollection;            IPixelBlock3 pixelblock3 = null;       

5、;      IRawBlocks  rawBlocks = IRawBlocksbands.Item0;            IRasterInfo  rasInfo = rawBlocks.RasterInfo;             long blockwidth = 0; &#

6、160;          long blockheight = 0;            System.Array pixels;            IPnt tlc = null;        &#

7、160;   object v;            double v1 = new double1000000;            int m = 0;            int lenght=0;   &#

8、160;        double, value=new double rasInfo.Width ,rasInfo.Height;            long bandCount = bands.Count;            do     

9、60;                      pixelblock3 = rasterCursor.PixelBlock as IPixelBlock3;                blockwidth = pixelblock3.Width

10、;                blockheight = pixelblock3.Height;                pixelblock3.Mask255;          

11、0;     /value = new doubleblockwidth, blockheight;                for int k = 0; k < bandCount; k+                  &

12、#160;                 /Get the pixel array.                    pixels = System.Arraypixelblock3.get_PixelDatak;  &#

13、160;                 for long i = 0; i < blockwidth; i+                            

14、;                for long j = 0; j < blockheight; j+                             &

15、#160;                      /Get the pixel value                        

16、0;       v = pixels.GetValuei, j;                            /Do something with the value./此处的数据提取是以64*64单元提取，需要主意顺序，该段代码前一讲中已说明。 

17、                           v1m =Convert .ToDouble v;                   

18、0;        lenght =m;                            m+;             &

19、#160;                                                  /Set the p

20、ixel array to the pixel block.                    pixelblock3.set_PixelDatak, pixels;                     

21、;           /Write back to the raster.                tlc = rasterCursor.TopLeft;                ras

22、terEdit.Writetlc, IPixelBlockpixelblock3;                        while rasterCursor.Next = true;            System.Runtime.InteropServi

23、ces.Marshal.ReleaseComObjectrasterEdit;                    m=0;/此处的循环为了将提取的数据进展重新排序            while m < lenght       

24、;                     for int i = 0; i < rasInfo.Height; i=i+64                    forint j=0;j<rasInfo

25、.Width ;j=j+64                        forint b=0;b<64;b+                     

26、       forint d=0;d<64;d+                                if i + d < rasInfo.Height && j + b < rasI

27、nfo.Width                                                 &#

28、160;                  valuej + b,i + d = v1m;                            &#

29、160;       m+;                                            &#

30、160;           return value;        2三角型面积计算，cellx，celly为栅格每一格所代表的实际长度public double AreaBasedouble cellX,double cellY,double cellZ1,double cellZ2,double cellZ3          &

31、#160;         double value1, value2, value3;             value1 = Math.SqrtcellX * cellX + cellZ1 - cellZ2 * cellZ1 - cellZ2;            value2 =

32、 Math.SqrtcellX * cellX + cellZ3 - cellZ2 * cellZ3 - cellZ2;            value3 = Math.SqrtcellX * cellX + cellY * cellY + cellZ3 - cellZ1 * cellZ3 - cellZ1;            double length=value1+v

33、alue2+value3;            return   Math.Sqrtlength * length -  value1 * length -  value2 * length -  value3;                   3总体面积计算：此

34、处采用一个sbyte三维数组，将每个栅格格子分成八个小三角形，目的是为了记录该快小三角形是否被之前相邻的格子使用过。以免造成面积重复计算public string  AreaCalculate                    try               &

35、#160;            IRasterBandCollection rbcollection = raster as IRasterBandCollection;                IRasterDataset rasterDataset1 = new RasterDatasetClass;  &

36、#160;             rasterDataset1 = rbcollection.Item0.RasterDataset as IRasterDataset;                IRasterDataset2 rasterDataset2 = rasterDataset1 as IRasterDatas

37、et2;                 double cellX, cellY;                int Xnum, Ynum;            &

38、#160;   RasterInfo rasterInfo = new RasterInfo;                rasterInfo.GetRasterInforaster;                Xnum = rasterInfo.Height; 

39、;               Ynum = rasterInfo.Width;                cellX = rasterInfo.xmax - rasterInfo.xmin / Xnum;        &#

40、160;       cellY = rasterInfo.ymax - rasterInfo.ymin / Ynum;                double, value = new doubleYnum, Xnum;             

41、0;  sbyte, , bool1 = new sbyteYnum, Xnum, 8;                for int i = 0; i < Ynum; i+                    for int j = 0;

42、j < Xnum; j+                        for int k = 0; k < 8; k+                    

43、;        bool1i, j, k = 0;                value = ReadWriteRawBlocksrasterDataset2;                double area = 0;

44、                double sarea = 0;                / return value0, 0;             

45、0;  for int i = 1; i < Ynum - 1; i+                                    for int j = 1; j < Xnum - 1; j+   

46、;                                         if valuei, j != -32768       

47、;                                             if valuei, j - 1 != -32768 && v

48、aluei + 1, j != -32768                                               &

49、#160;            if bool1i, j, 4 = 0 && bool1i, j, 5 = 0 && bool1i, j - 1, 7 = 0 && bool1i + 1, j, 2 = 0                   

50、60;                                                sarea = AreaBasecel

51、lX, cellY, valuei, j - 1, valuei, j, valuei + 1, j;                                    area = area + sarea;   

52、                                 bool1i, j, 4 = 1;               

53、                     bool1i, j, 5 = 1;                            

54、;        bool1i, j - 1, 7 = 1;                                    bool1i + 1, j, 2 =

55、1;                                                  

56、60;                                     if valuei + 1, j != -32768 && valuei + 1, j - 1 != -32768   &#

57、160;                                                  &

58、#160;     if bool1i, j, 5 = 0 && bool1i + 1, j, 2 = 0 && bool1i + 1, j, 3 = 0 && bool1i + 1, j - 1, 0 = 0                         &

59、#160;                                          sarea = AreaBasecellX, cellY, valuei, j, valuei + 1,

60、j, valuei + 1, j - 1;                                    area = area + sarea;        

61、                            bool1i, j, 5 = 1;                    

62、                bool1i + 1, j, 2 = 1;                                &

63、#160;   bool1i + 1, j, 3 = 1;                                    bool1i + 1, j - 1, 0 = 1;   &#

64、160;                                                  &

65、#160;                                 if valuei + 1, j != -32768 && valuei + 1, j + 1 != -32768       

66、                                                   

67、;  if bool1i, j, 6 = 0 && bool1i + 1, j, 0 = 0 && bool1i + 1, j, 1 = 0 && bool1i + 1, j + 1, 3 = 0                             

68、;                                       sarea = AreaBasecellX, cellY, valuei, j, valuei + 1, j, valuei + 1, j + 1;&

69、#160;                                   area = area + sarea;           

70、0;                        bool1i, j, 6 = 1;                       

71、0;            bool1i + 1, j, 0 = 1;                                   

72、 bool1i + 1, j, 1 = 1;                                    bool1i + 1, j + 1, 3 = 1;       

73、                                                   

74、;                              if valuei, j + 1 != -32768 && valuei + 1, j != -32768           &#

75、160;                                                if bool1i, j, 6 =

76、0 && bool1i, j, 7 = 0 && bool1i, j + 1, 4 = 0 && bool1i + 1, j, 1 = 0                                   &#

77、160;                                sarea = AreaBasecellX, cellY, valuei, j + 1, valuei, j, valuei + 1, j;      

78、60;                             area = area + sarea;                  

79、                  bool1i, j, 6 = 1;                              

80、      bool1i, j, 7 = 1;                                    bool1i, j + 1, 4 = 1;   &#

81、160;                                bool1i + 1, j, 1 = 1;              

82、0;                                                   &#

83、160;                                                   

84、                                                   

85、0;                                                double marea = cellX * ce

86、llY / 8;                for int i = 0; i < Ynum; i+                    for int j = 0; j < Xnum; j+    &#

87、160;                                       if valuei, j != -32768        &#

88、160;                                           for int k = 0; k < 8; k+   

89、60;                                                  &#

90、160;     if bool1i, j, k = 0                                    area = area + marea;   

91、60;                                                   &

92、#160;                                return area.ToString ;                

93、        catch Exception ex                           return ex.ToString ;           

94、          以下为主函数代码，写在mapcontrol的mousedown事件中4栅格裁剪：首先交互式绘制矢量多面性，采用IExtractionOp接口对原始栅格进展裁剪提取。IGeometry geometry = null;                geometry = axMapControl1.TrackPolygon; &#

95、160;              drawMapShape1geometry, 200, 50, 0;                IPolygon pPolygon = geometry as IPolygon;        

96、60;       IExtractionOp extraction = new RasterExtractionOpClass;  /此段代码为自己程序内部需要，自己可以直接将layer定义为要裁剪的栅格                ILayer layer = null;       

97、0;        for int i = 0; i < axMapControl1.LayerCount; i+                                   

98、 layer = axMapControl1.Map.get_Layeri as IRasterLayer;                    if this boBox1.SelectedItem.ToString = layer.Name             

99、0;                              break;                     &

100、#160;                  else                                &

101、#160;           MessageBox.Show"找¨°不?到Ì?对?应®|图ª?层?"                            &

102、#160;                          IRasterLayer rasterLayer = layer as IRasterLayer;                IRaster r

103、aster = rasterLayer.Raster;                           IRasterProps pProps = raster as IRasterProps;            

104、0;   object cellSizeProvider = pProps.MeanCellSize.X;                IGeoDataset pInputDataset = raster as IGeoDataset;                IExtr

105、actionOp pExtractionOp = new RasterExtractionOpClass;                IRasterAnalysisEnvironment pRasterAnaEnvir = pExtractionOp as IRasterAnalysisEnvironment;           &

106、#160;    pRasterAnaEnvir.SetCellSizeesriRasterEnvSettingEnum.esriRasterEnvValue, ref cellSizeProvider;                object extentProvider = pPolygon.Envelope;        &#

107、160;       object snapRasterData = Type.Missing;                pRasterAnaEnvir.SetExtentesriRasterEnvSettingEnum.esriRasterEnvValue, ref extentProvider, ref snapRasterData;   

108、0;            IGeoDataset pOutputDataset = pExtractionOp.PolygonpInputDataset, pPolygon, true;                IRaster newRaster = pOutputDataset as IRaster;  &#

109、160;              IRaster clipRaster;  /裁?切D后¨®得Ì?到Ì?的Ì?IRaster                if pOutputDataset is IRasterLayer   

110、                                 IRasterLayer rasterLayer1 = pOutputDataset as IRasterLayer;         

111、;           clipRaster = rasterLayer1.Raster;                                else if pOutputDataset

112、is IRasterDataset                                    IRasterDataset rasterDataset = pOutputDataset as IRasterDataset;  &

113、#160;                 clipRaster = rasterDataset.CreateDefaultRaster;                          

114、      else if pOutputDataset is IRaster                                    clipRaster = pOutputDataset as

115、 IRaster;                                else                  &#

116、160;                 return;                IRasterProps pProps = raster as IRasterProps;         &#

117、160;      object cellSizeProvider = pProps.MeanCellSize.X;                IGeoDataset pInputDataset = raster as IGeoDataset;                IExtractionOp pExtr