C# 使用 GDAL 读取 DEM.tif 文件并计算填方挖方

以下是使用 C# 代码读取 DEM.tif 文件并计算填方挖方的示例：

using System;
using System.IO;
using System.Linq;
using OSGeo.GDAL;

namespace FillCutVolume
{
    class Program
    {
        static void Main(string[] args)
        {
            // 注册 GDAL 驱动
            Gdal.AllRegister();

            // 打开 DEM 文件
            var demFile = 'dem.tif';
            var demDs = Gdal.Open(demFile, Access.GA_ReadOnly);
            if (demDs == null)
            {
                Console.WriteLine($'无法打开 DEM 文件 {demFile}');
                return;
            }

            // 读取 DEM 数据
            var cols = demDs.RasterXSize;
            var rows = demDs.RasterYSize;
            var band = demDs.GetRasterBand(1);
            var demData = new float[cols * rows];
            band.ReadRaster(0, 0, cols, rows, demData, cols, rows, 0, 0);

            // 计算填方挖方
            var fillVolume = 0.0;
            var cutVolume = 0.0;
            var cellSize = 1.0; // DEM 像元大小，单位为米
            var zFactor = 1.0; // DEM 高程缩放系数
            for (int i = 0; i < rows - 1; i++)
            {
                for (int j = 0; j < cols - 1; j++)
                {
                    var idx1 = i * cols + j;
                    var idx2 = (i + 1) * cols + j;
                    var idx3 = (i + 1) * cols + j + 1;
                    var idx4 = i * cols + j + 1;
                    var z1 = demData[idx1] * zFactor;
                    var z2 = demData[idx2] * zFactor;
                    var z3 = demData[idx3] * zFactor;
                    var z4 = demData[idx4] * zFactor;
                    var area = cellSize * cellSize;
                    var volume = (z1 + z2 + z3 + z4) / 4.0 * area;
                    if (z1 == 0 || z2 == 0 || z3 == 0 || z4 == 0)
                    {
                        // 如果四个点中有一个点高程为 0，则不计算填方挖方
                        continue;
                    }
                    var minZ = new[] { z1, z2, z3, z4 }.Min();
                    var maxZ = new[] { z1, z2, z3, z4 }.Max();
                    if (maxZ - minZ > 0)
                    {
                        // 如果四个点中的最大高程与最小高程不相等，则计算填方挖方
                        if (z1 == minZ)
                        {
                            fillVolume += volume;
                        }
                        else if (z1 == maxZ)
                        {
                            cutVolume += volume;
                        }
                        if (z2 == minZ)
                        {
                            fillVolume += volume;
                        }
                        else if (z2 == maxZ)
                        {
                            cutVolume += volume;
                        }
                        if (z3 == minZ)
                        {
                            fillVolume += volume;
                        }
                        else if (z3 == maxZ)
                        {
                            cutVolume += volume;
                        }
                        if (z4 == minZ)
                        {
                            fillVolume += volume;
                        }
                        else if (z4 == maxZ)
                        {
                            cutVolume += volume;
                        }
                    }
                }
            }

            // 输出填方挖方结果
            Console.WriteLine($'填方体积：{fillVolume} 立方米');
            Console.WriteLine($'挖方体积：{cutVolume} 立方米');

            // 关闭 DEM 文件
            demDs.Dispose();
            Gdal.GDALDestroyDriverManager();
        }
    }
}

其中，使用了 GDAL 库来读取 DEM 数据，计算填方挖方的方法是遍历 DEM 数据中的每个四边形，计算四边形的体积，并根据四个顶点的高程值判断该四边形是填方还是挖方。最终输出填方体积和挖方体积。需要注意的是，如果四个顶点中有一个顶点的高程值为 0，则不计算该四边形的填方挖方。