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

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

本文介绍使用GDAL库读取DEM.tif文件并计算填方挖方的方法，包括打开文件、获取文件信息、读取像素值、计算填方挖方等步骤，并提供完整代码示例。

1. 引入GDAL库头文件和命名空间

#include "gdal_priv.h"
using namespace std;

2. 打开dem.tif文件

GDALDataset* poDataset;
poDataset = (GDALDataset*)GDALOpen('dem.tif', GA_ReadOnly);
if (poDataset == NULL)
{
    cout << "Open failed." << endl;
exit(1);
}

3. 获取dem.tif文件的基本信息，如行数、列数、波段数等

int nXSize = poDataset->GetRasterXSize(); //列数
int nYSize = poDataset->GetRasterYSize(); //行数
int nBands = poDataset->GetRasterCount(); //波段数

4. 读取dem.tif文件的像素值

GDALRasterBand* poBand;
poBand = poDataset->GetRasterBand(1); //获取第1个波段
float* pafScanline = new float[nXSize]; //申请一行像素值的内存空间
for (int i = 0; i < nYSize; i++)
{
    poBand->RasterIO(GF_Read, 0, i, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0); //读取一行像素值
    //对一行像素值进行处理
}
delete[] pafScanline; //释放内存空间

5. 计算填方挖方

对于每个像素点，可以计算其高程值与基准高程值之差，若差值为正，则为填方，若差值为负，则为挖方。将所有像素点的填方和挖方相加即可得到总的填方和挖方。

float baseHeight = 0; //基准高程值
float totalFill = 0; //总的填方
float totalExcavate = 0; //总的挖方
for (int i = 0; i < nYSize; i++)
{
    poBand->RasterIO(GF_Read, 0, i, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0); //读取一行像素值
    for (int j = 0; j < nXSize; j++)
    {
        float height = pafScanline[j]; //当前像素点的高程值
        float diff = height - baseHeight; //当前像素点的高程值与基准高程值之差
        if (diff > 0)
            totalFill += diff; //填方
        else
            totalExcavate -= diff; //挖方
    }
}
cout << "Total fill: " << totalFill << endl;
cout << "Total excavate: " << totalExcavate << endl;

完整代码

#include <iostream>
#include "gdal_priv.h"
using namespace std;

int main()
{
    GDALAllRegister();

    GDALDataset* poDataset;
    poDataset = (GDALDataset*)GDALOpen('dem.tif', GA_ReadOnly);
if (poDataset == NULL)
    {
        cout << "Open failed." << endl;
exit(1);
    }

    int nXSize = poDataset->GetRasterXSize(); //列数
    int nYSize = poDataset->GetRasterYSize(); //行数
    int nBands = poDataset->GetRasterCount(); //波段数

    GDALRasterBand* poBand;
    poBand = poDataset->GetRasterBand(1); //获取第1个波段
    float* pafScanline = new float[nXSize]; //申请一行像素值的内存空间

    float baseHeight = 0; //基准高程值
    float totalFill = 0; //总的填方
    float totalExcavate = 0; //总的挖方
    for (int i = 0; i < nYSize; i++)
    {
        poBand->RasterIO(GF_Read, 0, i, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0); //读取一行像素值
        for (int j = 0; j < nXSize; j++)
        {
            float height = pafScanline[j]; //当前像素点的高程值
            float diff = height - baseHeight; //当前像素点的高程值与基准高程值之差
            if (diff > 0)
                totalFill += diff; //填方
            else
                totalExcavate -= diff; //挖方
        }
    }
    cout << "Total fill: " << totalFill << endl;
    cout << "Total excavate: " << totalExcavate << endl;

delete[] pafScanline;
    GDALClose(poDataset);
    return 0;
}