ARCGIS PRO DSK 利用两期地表DEM数据计算工程土方量

利用两期地表DEM数据计算工程土方量需要准许以下数据：
当前地图有3个图层，两个栅格图层和一个矢量图层
两个栅格图层：beforeDem为工程施工前的地表DEM模型
                         afterDem为工程施工后的地表DEM模型
一个矢量图层：pFeatureLayer
第一步生成掩膜：面转栅格 conversion.PolygonToRaster 生产名称为“poly”的栅格数据。
                             OID为pFeatureLayer的字段

 pFeatureLayer = PMap.FindLayers(ComboBox4.Text).FirstOrDefaultIf pFeatureLayer IsNot Nothing ThenGoTo qqwwDim maskRaster As String = Path.Combine(Project.Current.DefaultGeodatabasePath, $"poly")va = Geoprocessing.MakeValueArray(pFeatureLayer, "OID", maskRaster)  ' 直接传递环境参数列表Dim result = Await Geoprocessing.ExecuteToolAsync("PolygonToRaster_conversion", va)If result.IsFailed ThenThrow New Exception($"掩膜生成失败: {String.Join(vbCrLf, result.ErrorMessages)}")End If

第二步生成掩膜：按掩膜提取分别生成两个栅格图层
                          before_clipped：为工程施工前的地表DEM的提取模型
                          after_clipped   ：为工程施工后的地表DEM的提取模型

 Dim outputName As StringDim inputRaster As StringDim outputPath As StringDim envSettings = New List(Of KeyValuePair(Of String, Object))For i = 1 To 2If i = 1 Then inputRaster = beforeDem : outputName = "before_clipped"If i = 2 Then inputRaster = afterDem : outputName = "after_clipped"'直接传递环境参数列表envSettings = New List(Of KeyValuePair(Of String, Object)) From {New KeyValuePair(Of String, Object)("cellSize", "MAXOF"), New KeyValuePair(Of String, Object)("template", inputRaster)}outputPath = Path.Combine(Project.Current.DefaultGeodatabasePath, outputName)Dim extractParams = Geoprocessing.MakeValueArray(inputRaster, maskRaster, outputPath, "INSIDE")result = Await Geoprocessing.ExecuteToolAsync("ExtractByMask", extractParams)                      'ExtractByMaskIf result.IsFailed ThenMsgBox(inputRaster & "栅格按掩膜提取失败")End IfNext

第三步CutFill分析：执行CutFill分析，分析结果生成分析栅格数据CutFill_Result

            '1. 裁剪DEM到指定范围Dim clippedBefore = Path.Combine(Project.Current.DefaultGeodatabasePath, "before_clipped")Dim clippedAfter = Path.Combine(Project.Current.DefaultGeodatabasePath, "after_clipped")'2. 执行CutFill分析outputName = "CutFill_Result"outputPath = System.IO.Path.Combine(Project.Current.DefaultGeodatabasePath, outputName)Dim parameters = Geoprocessing.MakeValueArray(clippedBefore, clippedAfter, outputPath, "METER")result = Await Geoprocessing.ExecuteToolAsync("CutFill", parameters)If result.IsFailed ThenThrow New Exception($"CutFill分析失败: {String.Join(vbCrLf, result.ErrorMessages)}")End If

第四步分析成果按属性提取：ExtractByAttributes，提取的数据仍然为栅格数据
                                               属性提取Where 子句使用 SQL 查询 
                                               "VOLUME < 0"  填方
                                               "VOLUME > 0"  挖方
                                               "VOLUME = 0"  未变化
       输出栅格的属性表用于确定显示方式，并且分别将正体积和负体积视为挖出材料的位置（已移除）和填充材料的位置（已添加）。

' 计算挖方量
Dim cutParams = Geoprocessing.MakeValueArray(cutFillRasterPath, "VOLUME > 0", GeodaPath & "\WFMJ")    '按属性提取
gpResult = Await Geoprocessing.ExecuteToolAsync("ExtractByAttributes", cutParams)
If gpResult.IsFailed ThenMsgBox("挖方区域提取失败.")
End If' 计算填方量Dim fillParams = Geoprocessing.MakeValueArray(cutFillRasterPath, "VOLUME < 0", GeodaPath & "\" & "TFMJ")   '按属性提取gpResult = Await Geoprocessing.ExecuteToolAsync("ExtractByAttributes", fillParams)If gpResult.IsFailed ThenMsgBox("填方区域提取失败.")End If' 计算未变化方量Dim wbhParams = Geoprocessing.MakeValueArray(GeodaPath & "\" & cutFillRasterPath, "VOLUME = 0", GeodaPath & "\" & "WBHMJ")gpResult = Await Geoprocessing.ExecuteToolAsync("ExtractByAttributes_management", wbhParams)If gpResult.IsFailed ThenMsgBox("未变化区域提取失败.")End If

第五步获取栅格数据属性如：

• x 方向上的像元大小—x 方向上的像元大小。
• y 方向上的像元大小—y 方向上的像元大小。

 va = Geoprocessing.MakeValueArray("TFMJ", {"CELLSIZEX", "CELLSIZEY"})gpResult = Await Geoprocessing.ExecuteToolAsync("management.GetRasterProperties", va)If gpResult.IsFailed ThenMsgBox("栅格数据属性提取失败.")End If

第六步分析栅格数据属性表访问Attributes：

 TryAwait QueuedTask.Run(Sub()' 获取栅格数据集Dim raster As Raster = rasterLayer.GetRaster()If raster Is Nothing Then Return' 正确检查属性表的方法Dim hasTable As Boolean = FalseTry' 尝试获取属性表Using table As Table = raster.GetAttributeTable()hasTable = True' 获取字段信息'Dim fields As IReadOnlyList(Of Field) = table.GetDefinition().GetFields()'Dim fieldNames As String = String.Join(", ", fields.Select(Function(f) f.Name))'MessageBox.Show($"找到属性表，包含字段: {fieldNames}")Dim tabrow = table.Search(Nothing, False)Do While tabrow.MoveNext'访问属性表 LoopEnd UsingCatch ex As Exception When TypeOf ex Is NotSupportedException OrElseTypeOf ex Is InvalidOperationException' 捕获不支持属性表的异常hasTable = FalseEnd TryIf Not hasTable ThenMessageBox.Show("该栅格没有属性表")End IfEnd Sub)Catch ex As ExceptionMessageBox.Show($"访问属性表时出错: {ex.Message}")End Try

如果只需要计算方量，第四步、第五步可省略。

运行后，可得：