最新AE二次开发进行DEM表面积求算.docx

1、最新AE二次开发进行DEM表面积求算AE二次开发进行DEM表面积求算四、算法代码此后为自定义类中的代码（1）栅格数据提取：采用二维数组对栅格数据进行结果存储，代码为arcengine帮助文档内部搜索出然后修正的public double, ReadWriteRawBlocks(IRasterDataset2 rasterDs) /Create a raster. IRaster2 raster2 = rasterDs.CreateFullRaster() as IRaster2; /Create a raster cursor with a system-optimized pixel blo

2、ck size by passing a null. IRasterCursor rasterCursor = raster2.CreateCursorEx(null); /Use the IRasterEdit interface. IRasterEdit rasterEdit = raster2 as IRasterEdit; /Loop through each band and pixel block. IRasterBandCollection bands = rasterDs as IRasterBandCollection; IPixelBlock3 pixelblock3 =

3、null; IRawBlocks rawBlocks = (IRawBlocks)bands.Item(0); IRasterInfo rasInfo = rawBlocks.RasterInfo; long blockwidth = 0; long blockheight = 0; System.Array pixels; IPnt tlc = null; object v; double v1 = new double1000000; int m = 0; int lenght=0; double, value=new double rasInfo.Width ,rasInfo.Heigh

4、t; long bandCount = bands.Count; do pixelblock3 = rasterCursor.PixelBlock as IPixelBlock3; blockwidth = pixelblock3.Width; blockheight = pixelblock3.Height; pixelblock3.Mask(255); /value = new doubleblockwidth, blockheight; for (int k = 0; k bandCount; k+) /Get the pixel array. pixels = (System.Arra

5、y)pixelblock3.get_PixelData(k); for (long i = 0; i blockwidth; i+) for (long j = 0; j blockheight; j+) /Get the pixel value v = pixels.GetValue(i, j); /Do something with the value./此处的数据提取是以64*64单元提取，需要主意顺序，该段代码前一讲中已说明。 v1m =Convert .ToDouble (v); lenght =m; m+; /Set the pixel array to the pixel blo

6、ck. pixelblock3.set_PixelData(k, pixels); /Write back to the raster. tlc = rasterCursor.TopLeft; rasterEdit.Write(tlc, (IPixelBlock)pixelblock3); while (rasterCursor.Next() = true); System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterEdit); m=0;/此处的循环为了将提取的数据进行重新排序 while (m lenght) for (in

7、t i = 0; i rasInfo.Height; i=i+64) for(int j=0;jrasInfo .Width ;j=j+64) for(int b=0;b64;b+) for(int d=0;d64;d+) if (i + d rasInfo.Height & j + b rasInfo.Width) valuej + b,i + d = v1m; m+; return value; （2）三角型面积计算，cellx，celly为栅格每一格所代表的实际长度public double AreaBase(double cellX,double cellY,double cellZ1

8、,double cellZ2,double cellZ3) double value1, value2, value3; value1 = Math.Sqrt(cellX * cellX + (cellZ1 - cellZ2) * (cellZ1 - cellZ2); value2 = Math.Sqrt(cellX * cellX + (cellZ3 - cellZ2) * (cellZ3 - cellZ2); value3 = Math.Sqrt(cellX * cellX + cellY * cellY + (cellZ3 - cellZ1) * (cellZ3 - cellZ1); d

9、ouble length=value1+value2+value3; return Math.Sqrt(length * (length - value1) * (length - value2) * (length - value3); （3）总体面积计算：此处采用一个sbyte三维数组，将每个栅格格子分成八个小三角形，目的是为了记录该快小三角形是否被之前相邻的格子使用过。以免造成面积重复计算public string AreaCalculate() try IRasterBandCollection rbcollection = raster as IRasterBandCollectio

10、n; IRasterDataset rasterDataset1 = new RasterDatasetClass(); rasterDataset1 = rbcollection.Item(0).RasterDataset as IRasterDataset; IRasterDataset2 rasterDataset2 = rasterDataset1 as IRasterDataset2; double cellX, cellY; int Xnum, Ynum; RasterInfo rasterInfo = new RasterInfo(); rasterInfo.GetRasterI

11、nfo(raster); Xnum = rasterInfo.Height; Ynum = rasterInfo.Width; cellX = (rasterInfo.xmax - rasterInfo.xmin) / Xnum; cellY = (rasterInfo.ymax - rasterInfo.ymin) / Ynum; double, value = new doubleYnum, Xnum; sbyte, , bool1 = new sbyteYnum, Xnum, 8; for (int i = 0; i Ynum; i+) for (int j = 0; j Xnum; j

12、+) for (int k = 0; k 8; k+) bool1i, j, k = 0; value = ReadWriteRawBlocks(rasterDataset2); double area = 0; double sarea = 0; / return value0, 0; for (int i = 1; i Ynum - 1; i+) for (int j = 1; j Xnum - 1; j+) if (valuei, j != -32768) if (valuei, j - 1 != -32768 & valuei + 1, j != -32768) if (bool1i,

13、 j, 4 = 0 & bool1i, j, 5 = 0 & bool1i, j - 1, 7 = 0 & bool1i + 1, j, 2 = 0) sarea = AreaBase(cellX, cellY, valuei, j - 1, valuei, j, valuei + 1, j); area = area + sarea; bool1i, j, 4 = 1; bool1i, j, 5 = 1; bool1i, j - 1, 7 = 1; bool1i + 1, j, 2 = 1; if (valuei + 1, j != -32768 & valuei + 1, j - 1 !=

14、 -32768) if (bool1i, j, 5 = 0 & bool1i + 1, j, 2 = 0 & bool1i + 1, j, 3 = 0 & bool1i + 1, j - 1, 0 = 0) sarea = AreaBase(cellX, cellY, valuei, j, valuei + 1, j, valuei + 1, j - 1); area = area + sarea; bool1i, j, 5 = 1; bool1i + 1, j, 2 = 1; bool1i + 1, j, 3 = 1; bool1i + 1, j - 1, 0 = 1; if (valuei

15、 + 1, j != -32768 & valuei + 1, j + 1 != -32768) if (bool1i, j, 6 = 0 & bool1i + 1, j, 0 = 0 & bool1i + 1, j, 1 = 0 & bool1i + 1, j + 1, 3 = 0) sarea = AreaBase(cellX, cellY, valuei, j, valuei + 1, j, valuei + 1, j + 1); area = area + sarea; bool1i, j, 6 = 1; bool1i + 1, j, 0 = 1; bool1i + 1, j, 1 =

16、 1; bool1i + 1, j + 1, 3 = 1; if (valuei, j + 1 != -32768 & valuei + 1, j != -32768) if (bool1i, j, 6 = 0 & bool1i, j, 7 = 0 & bool1i, j + 1, 4 = 0 & bool1i + 1, j, 1 = 0) sarea = AreaBase(cellX, cellY, valuei, j + 1, valuei, j, valuei + 1, j); area = area + sarea; bool1i, j, 6 = 1; bool1i, j, 7 = 1

17、; bool1i, j + 1, 4 = 1; bool1i + 1, j, 1 = 1; double marea = cellX * cellY / 8; for (int i = 0; i Ynum; i+) for (int j = 0; j Xnum; j+) if (valuei, j != -32768) for (int k = 0; k 8; k+) if (bool1i, j, k = 0) area = area + marea; return area.ToString (); catch (Exception ex) return ex.ToString (); 以下为主函数代码，写在mapcontrol的mousedown事件中（4）栅格裁剪：首先交互式绘制矢量多面性，采用IExtractionOp接口对原始栅格进行裁剪提取。IGeometry geometry = null; geometry = axMapControl1.TrackPolygon(); drawMapShape1(geometry, 200, 50, 0); IPolygon pPolygon = geometry as IPolygon; IExtractionOp extraction = new RasterExtractionOpClass();