前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >GDAL入门-使用GDAL进行遥感影像NDVI的计算(C++版本)

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

作者头像
卡尔曼和玻尔兹曼谁曼
发布2019-01-22 10:58:07
4.1K1
发布2019-01-22 10:58:07
举报
文章被收录于专栏:给永远比拿愉快

不知道该怎么详细地写,直接上代码算了! (开发环境的搭建参考我的博文:GDAL开发环境搭建(VS2010 C++版)

代码语言:javascript
复制
#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;
}
本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2016年09月07日,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档