C语言实现读取DEM数据并计算填方挖方

由于DEM数据是二维数组形式的高程数据，因此我们可以通过读取DEM数据的方式来计算填方挖方。

具体实现步骤如下：

1. 读取DEM数据，可以使用GDAL库来读取各种格式的DEM数据。

2. 计算DEM数据的最小高程和最大高程，即DEM数据的高程范围。

3. 根据高程范围，将高程值划分为若干个等高线段，每个等高线段的高程范围为一个固定的值。这里可以根据实际情况来确定等高线段的高程范围，比如每10米一个等高线段。

4. 对于每个等高线段，计算该等高线段下方的面积和上方的面积，即填方和挖方的面积。

5. 将每个等高线段的填方和挖方面积累加起来，即可得到总的填方和挖方面积。

6. 根据填方和挖方面积的差值，即可计算出填方和挖方的体积。

下面是一份简单的代码示例，演示了如何读取DEM数据并计算填方挖方：

#include <stdio.h>
#include <gdal.h>

int main()
{
    // 打开DEM数据文件
    GDALDatasetH hDataset = GDALOpen('dem.tif', GA_ReadOnly);
    if (hDataset == NULL)
    {
        printf('Open DEM file failed!\n');
        return -1;
    }

    // 获取DEM数据的基本信息
    int nXSize = GDALGetRasterXSize(hDataset);
    int nYSize = GDALGetRasterYSize(hDataset);
    double adfGeoTransform[6];
    GDALGetGeoTransform(hDataset, adfGeoTransform);

    // 读取DEM数据
    float *pafData = new float[nXSize * nYSize];
    GDALRasterBandH hBand = GDALGetRasterBand(hDataset, 1);
    GDALRasterIO(hBand, GF_Read, 0, 0, nXSize, nYSize, pafData, nXSize, nYSize, GDT_Float32, 0, 0);

    // 计算DEM数据的高程范围
    double fMinHeight = pafData[0];
    double fMaxHeight = pafData[0];
    for (int i = 0; i < nXSize * nYSize; i++)
    {
        if (pafData[i] < fMinHeight)
            fMinHeight = pafData[i];
        if (pafData[i] > fMaxHeight)
            fMaxHeight = pafData[i];
    }

    // 划分等高线段
    double fContourInterval = 10.0; // 每10米一个等高线段
    int nContourNum = (int)((fMaxHeight - fMinHeight) / fContourInterval) + 1;
    double *pafContourHeight = new double[nContourNum];
    for (int i = 0; i < nContourNum; i++)
    {
        pafContourHeight[i] = fMinHeight + fContourInterval * i;
    }

    // 计算填方挖方
    double fFillVolume = 0.0;
    double fExcavateVolume = 0.0;
    for (int i = 0; i < nContourNum - 1; i++)
    {
        double fContourHeight1 = pafContourHeight[i];
        double fContourHeight2 = pafContourHeight[i + 1];
        double fAreaAbove = 0.0;
        double fAreaBelow = 0.0;
        for (int j = 0; j < nXSize * nYSize; j++)
        {
            if (pafData[j] >= fContourHeight1 && pafData[j] < fContourHeight2)
            {
                double x = adfGeoTransform[0] + adfGeoTransform[1] * (j % nXSize) + adfGeoTransform[2] * (j / nXSize);
                double y = adfGeoTransform[3] + adfGeoTransform[4] * (j % nXSize) + adfGeoTransform[5] * (j / nXSize);
                if (pafData[j] >= (fContourHeight1 + fContourHeight2) / 2.0)
                    fAreaAbove += adfGeoTransform[1] * adfGeoTransform[5];
                else
                    fAreaBelow += adfGeoTransform[1] * adfGeoTransform[5];
            }
        }
        fFillVolume += fAreaAbove * fContourInterval;
        fExcavateVolume += fAreaBelow * fContourInterval;
    }

    // 输出填方挖方结果
    printf('Fill volume: %.2lf cubic meters\n', fFillVolume);
    printf('Excavate volume: %.2lf cubic meters\n', fExcavateVolume);

    // 释放内存
    delete[] pafData;
    delete[] pafContourHeight;
    GDALClose(hDataset);

    return 0;
}