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

本是古典何须时尚 2022-07-18 06:19 511阅读 0赞

不知道该怎么详细地写，直接上代码算了！  
（开发环境的搭建参考我的博文：[GDAL开发环境搭建（VS2010 C++版）][GDAL_VS2010 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;
    }

[GDAL_VS2010 C]: http://blog.csdn.net/theonegis/article/details/25590269