GDAL入门-使用GDAL进行遥感影像NDVI的计算（C++版本）

#include <iostream>
#include <array>
#include "gdal_priv.h"
#include "ogr_spatialref.h"


using std::cout;
using std::endl;
using std::array;

/*
    @beief 计算NDVI NDVI=(Red-NIR)/(Red+NIR)=(Band3-Band4)/(Band3+Band4)
    @param inputFileNames 输入参数，红波段和红外波段的两幅遥感影像的全路径。 如果影像是合成影像，则程序需要做相应的修改。
    @param outputFileName 输出参数，输出NDVI影像的全路径
*/
int NDVI(array<const char*, 2>& inputFileNames, const char* outputFileName)
{
    array<GDALDataset*, 2> inputDatasets;
    GDALDataset* outputDataset;

    GDALAllRegister();

    cout << "正在读取影像..." << endl;
    for (int i = 0; i < 2; i++)
    {
        inputDatasets[i] = (GDALDataset *)GDALOpen(inputFileNames[i], GA_ReadOnly);
        if (inputDatasets[i] == NULL)
        {
            cout << "影像读取失败：（" << inputFileNames[i] << "）！" << '\n';
            return EXIT_FAILURE;
        }
    }


    int imgSizeX = inputDatasets[0]->GetRasterXSize(); // 影像的宽度（像元数目）
    int imgSizeY = inputDatasets[0]->GetRasterYSize(); // 影像的高度（像元数目）

    const char* imageFormat = "GTiff";
    GDALDriver* gdalDriver = GetGDALDriverManager()->GetDriverByName(imageFormat);
    if (gdalDriver == NULL)
    { 
        cout << "创建输出影像失败！" << '\n';
        return EXIT_FAILURE;
    }

    outputDataset = gdalDriver->Create(outputFileName, imgSizeX, imgSizeY, 1, GDT_Float32, NULL);

    // 获取输入数据的地理变化信息
    double goeInformation[6];
    inputDatasets[0]->GetGeoTransform(goeInformation);
    // 读取输入数据的地理信息并写入输出文件
    const char* gdalProjection = inputDatasets[0]->GetProjectionRef();
    outputDataset->SetGeoTransform(goeInformation);
    outputDataset->SetProjection(gdalProjection);

    cout << "正在进行数据处理..." << '\n';
    // 取得红波段和近红外波段
    // 读取数据，因为只有一个波段，所以Get第一波段的数据。如果是合成影像，这里依次读取各个波段的值（计数从1开始）。
    GDALRasterBand* raseterBandRed = inputDatasets[0]->GetRasterBand(1);
    GDALRasterBand* raseterBandNIR = inputDatasets[1]->GetRasterBand(1);
    GDALRasterBand* outputRasterBand = outputDataset->GetRasterBand(1);
    // 申请存储空间，为一行的大小
    float* bufferBlockRed = (float*)CPLMalloc(sizeof(float) * imgSizeX);
    float* bufferBlockNIR = (float*)CPLMalloc(sizeof(float) * imgSizeX);
    float* outputBufferBlock = (float*)CPLMalloc(sizeof(float) * imgSizeX);

    // 进行NDWI的计算
    for (int i = 0; i < imgSizeY; i++)
    {
        raseterBandRed->RasterIO(GF_Read, 0, i, imgSizeX, 1, bufferBlockRed, imgSizeX, 1, GDT_UInt32, 0, 0);
        raseterBandNIR->RasterIO(GF_Read, 0, i, imgSizeX, 1, bufferBlockNIR, imgSizeX, 1, GDT_UInt32, 0, 0);
        for (int j = 0; j < imgSizeX; j++)
            outputBufferBlock[j] = (bufferBlockRed[j] - bufferBlockNIR[j]) / (bufferBlockRed[j] + bufferBlockNIR[j]);
        outputRasterBand->RasterIO(GF_Write, 0, i, imgSizeX, 1, outputBufferBlock, imgSizeX, 1, GDT_Float32, 0, 0); // 写入数据
    }

    // 释放资源
    CPLFree(bufferBlockRed);
    CPLFree(bufferBlockNIR);
    CPLFree(outputBufferBlock);
    GDALClose(inputDatasets[0]);
    GDALClose(inputDatasets[1]);
    GDALClose(outputDataset);

    cout << "处理结束！" << endl;

    return EXIT_SUCCESS;
}

int main()
{
    array<const char*, 2> inputs = {
        "C:/Users/theone/Desktop/LT51240392010131BKT01_B3.TIF",
        "C:/Users/theone/Desktop/LT51240392010131BKT01_B4.TIF" };
    const char* outputs = "C:/Users/theone/Desktop/LT51240392010131BKT01_NDVI.TIF";
    int isSuccess = NDVI(inputs, outputs);
    getchar(); // 防止控制台一闪而过
    return isSuccess;
}