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

本文链接：https://blog.csdn.net/xa13850869/article/details/147189527

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

 pFeatureLayer = PMap.FindLayers(ComboBox4.Text).FirstOrDefault
 If pFeatureLayer IsNot Nothing Then
     GoTo qqww
     Dim 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 Then
         Throw New Exception($"掩膜生成失败: {String.Join(vbCrLf, result.ErrorMessages)}")
     End If

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

 Dim outputName As String
 Dim inputRaster As String
 Dim outputPath As String
 Dim envSettings = New List(Of KeyValuePair(Of String, Object))
 For i = 1 To 2
     If 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)                      'ExtractByMask
     If result.IsFailed Then
         MsgBox(inputRaster & "栅格按掩膜提取失败")
     End If
 Next

第三步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 Then
                Throw 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 Then
    MsgBox("挖方区域提取失败.")
End If
 ' 计算填方量
 Dim fillParams = Geoprocessing.MakeValueArray(cutFillRasterPath, "VOLUME < 0", GeodaPath & "\" & "TFMJ")   '按属性提取
 gpResult = Await Geoprocessing.ExecuteToolAsync("ExtractByAttributes", fillParams)
 If gpResult.IsFailed Then
    MsgBox("填方区域提取失败.")
 End If
 ' 计算未变化方量
 Dim wbhParams = Geoprocessing.MakeValueArray(GeodaPath & "\" & cutFillRasterPath, "VOLUME = 0", GeodaPath & "\" & "WBHMJ")
 gpResult = Await Geoprocessing.ExecuteToolAsync("ExtractByAttributes_management", wbhParams)
 If gpResult.IsFailed Then
     MsgBox("未变化区域提取失败.")
 End If

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

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

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

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

 Try
     Await QueuedTask.Run(Sub()
                              ' 获取栅格数据集
                              Dim raster As Raster = rasterLayer.GetRaster()
                              If raster Is Nothing Then Return
                              ' 正确检查属性表的方法
                              Dim hasTable As Boolean = False
                              Try
                                  ' 尝试获取属性表
                                  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
                                         '访问属性表 
                                      Loop
                                  End Using
                              Catch ex As Exception When TypeOf ex Is NotSupportedException OrElse
                                TypeOf ex Is InvalidOperationException
                                  ' 捕获不支持属性表的异常
                                  hasTable = False
                              End Try

                              If Not hasTable Then
                                  MessageBox.Show("该栅格没有属性表")
                              End If
                          End Sub)
 Catch ex As Exception
     MessageBox.Show($"访问属性表时出错: {ex.Message}")
 End Try

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

运行后，可得：