AE二次開發(fā)進行DEM表面積求算

1、.四、算法代碼此后為自定義類中的代碼（1）柵格數(shù)據(jù)提?。翰捎枚S數(shù)組對柵格數(shù)據(jù)進行結(jié)果存儲，代碼為arcengine幫助文檔內(nèi)部搜索出然后修正的public double, ReadWriteRawBlocks(IRasterDataset2 rasterDs) /Create a raster. IRaster2 raster2 = rasterDs.CreateFullRaster() as IRaster2; /Create a raster cursor with a system-optimized pixel block size by passing a null. IRaste

2、rCursor 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 = null; IRawBlocks rawBlocks = (IRa

3、wBlocks)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.Height; long bandCount = bands.Count;

4、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.Array)pixelblock3.get_PixelData(k); f

5、or (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./此處的數(shù)據(jù)提取是以64*64單元提取，需要主意順序，該段代碼前一講中已說明。 v1m =Convert .ToDouble (v); lenght =m; m+; /Set the pixel array to the pixel block. pixelblock3.set_PixelData(k,

6、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;/此處的循環(huán)為了將提取的數(shù)據(jù)進行重新排序 while (m lenght) for (int i = 0; i rasInfo.Height; i=i+64

7、) 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,double cellZ2,double cellZ3) dou

8、ble 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); double length=value1+value2+value3

9、; return Math.Sqrt(length * (length - value1) * (length - value2) * (length - value3); （3）總體面積計算：此處采用一個sbyte三維數(shù)組，將每個柵格格子分成八個小三角形，目的是為了記錄該快小三角形是否被之前相鄰的格子使用過。以免造成面積重復(fù)計算public string AreaCalculate() try IRasterBandCollection rbcollection = raster as IRasterBandCollection; IRasterDataset rasterDataset1

10、= 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.GetRasterInfo(raster); Xnum = rasterInfo.He

11、ight; 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+) for (int k = 0; k 8; k+) bool1

12、i, 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, j, 4 = 0 & bool1i, j, 5 = 0 & bo

13、ol1i, 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 != -32768) if (bool1i, j, 5 = 0 & b

14、ool1i + 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 + 1, j != -32768 & valuei + 1, j

15、 + 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 = 1; bool1i + 1, j + 1, 3 = 1; if

16、(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; bool1i, j + 1, 4 = 1; bool1i +

17、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 (); 以下為主函數(shù)代碼，寫在mapcontrol的mousedown事件中（

18、4）柵格裁剪：首先交互式繪制矢量多面性，采用IExtractionOp接口對原始柵格進行裁剪提取。IGeometry geometry = null; geometry = axMapControl1.TrackPolygon(); drawMapShape1(geometry, 200, 50, 0); IPolygon pPolygon = geometry as IPolygon; IExtractionOp extraction = new RasterExtractionOpClass();/此段代碼為自己程序內(nèi)部需要，自己可以直接將layer定義為要裁剪的柵格 ILayer lay

19、er = null; for (int i = 0; i axMapControl1.LayerCount; i+) layer = axMapControl1.Map.get_Layer(i) as IRasterLayer; if (boBox1.SelectedItem.ToString() = layer.Name) break; else MessageBox.Show(找不?到?對?應(yīng)|圖?層?); IRasterLayer rasterLayer = layer as IRasterLayer; IRaster raster = rasterLayer.Raster; IRast

20、erProps pProps = raster as IRasterProps; object cellSizeProvider = pProps.MeanCellSize().X; IGeoDataset pInputDataset = raster as IGeoDataset; IExtractionOp pExtractionOp = new RasterExtractionOpClass(); IRasterAnalysisEnvironment pRasterAnaEnvir = pExtractionOp as IRasterAnalysisEnvironment; pRaste

21、rAnaEnvir.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref cellSizeProvider); object extentProvider = pPolygon.Envelope; object snapRasterData = Type.Missing; pRasterAnaEnvir.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, ref extentProvider, ref snapRasterData); IGeoDataset pOutp

22、utDataset = pExtractionOp.Polygon(pInputDataset, pPolygon, true); IRaster newRaster = pOutputDataset as IRaster; IRaster clipRaster; /裁?切D后得?到?的?IRaster if (pOutputDataset is IRasterLayer) IRasterLayer rasterLayer1 = pOutputDataset as IRasterLayer; clipRaster = rasterLayer1.Raster; else if (pOutputD

23、ataset is IRasterDataset) IRasterDataset rasterDataset = pOutputDataset as IRasterDataset; clipRaster = rasterDataset.CreateDefaultRaster(); else if (pOutputDataset is IRaster) clipRaster = pOutputDataset as IRaster; else return; IRasterProps pProps = raster as IRasterProps; object cellSizeProvider

24、= pProps.MeanCellSize().X; IGeoDataset pInputDataset = raster as IGeoDataset; IExtractionOp pExtractionOp = new RasterExtractionOpClass(); IRasterAnalysisEnvironment pRasterAnaEnvir = pExtractionOp as IRasterAnalysisEnvironment; pRasterAnaEnvir.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref cellSizeProvider); object extentProvider = pPolygon.Envelope; object snapRasterData = Type.Mi