又经过了几天的学习，完成了第二个实验，相较于第一个实验，第二个实验在进行的过程中确实遇到了一些问题，耗时也更长，不过总归是完成了。

任务

计算植被指数，并通过阈值法提取植被

代码

注：部分函数的注释由于在第一篇帖子里出现了，这里就不赘言了
程序的结构如下：

主函数proj2.cpp

/*
* 用于学习GDAL编程 实验二 计算植被指数，并通过阈值法提取植被
* 时间 20210705
* 本次实验采用的影像为四波段的影像
*/#include
#include "gdal.h"
#include "gdal_priv.h"
#include "Func.h"using namespace std;int main()
{// 读取影像char* openpath = "D:\\Practice\\org\\exp2.tif";GDALDataset* mDataset;GDALAllRegister();// 只读方式打开mDataset = (GDALDataset*)GDALOpen(openpath, GA_ReadOnly);if (mDataset == NULL){cout << "未能正确加载数据!!!" << endl;system("pause");GDALDestroyDriverManager();}// 获取影像数据int mX = mDataset->GetRasterXSize();int mY = mDataset->GetRasterYSize();GDALDataType mImgDataType = mDataset->GetRasterBand(1)->GetRasterDataType();// 获取影像的地理信息和投影信息double geoInformation[6];mDataset->GetGeoTransform(geoInformation);const char* gdalProjection = mDataset->GetProjectionRef();// 新建一个驱动，影像格式为GTiffGDALDriver* hDriver = GetGDALDriverManager()->GetDriverByName("GTiff");// 调用CalcuNDVI 进行计算，需要先建立一个数据集保存NDVI结果,NDVI中的数据类型是floatchar* savepath = "D:\\Practice\\res\\exp2_ndvi.tif";GDALDataset* mNDVIset = hDriver->Create(savepath, mX, mY, 1, GDT_Float32, NULL);cout << "正在计算NDVI" << endl;CalcuNDVI(mDataset, mNDVIset);// 将影像的地理信息写到NDVI数据集中mNDVIset->SetGeoTransform(geoInformation);mNDVIset->SetProjection(gdalProjection);// 调用ExtractVegetation进行提取，需要先建立一个数据集保存提取结果char* savepath2 = "D:\\Practice\\res\\exp2_vegetation.tif";GDALDataset* mVegetationSet = hDriver->Create(savepath2, mX, mY, 4, mImgDataType, NULL);cout << "正在提取植被区域！！！" << endl;ExtractVegetation(mDataset, mNDVIset, mVegetationSet, 0.35, 0.7);// 将影像的地理信息写到植被数据集中mVegetationSet->SetGeoTransform(geoInformation);mVegetationSet->SetProjection(gdalProjection);// 关闭数据集，释放内存GDALClose(mDataset);GDALClose(mNDVIset);GDALClose(mVegetationSet);GDALDestroyDriverManager();cout << "处理结束！！！" << endl;getchar();return 0;
}

头文件Func.h

#pragma once
#include 
#include "gdal.h"
#include "gdal_priv.h"
#include "ogr_spatialref.h"using namespace std;// 计算NDVI,mDataset为输入数据，mNDVI为计算结果
void CalcuNDVI(GDALDataset* mDataset,GDALDataset* &mNDVI);// 利用NDVI，采用阈值分割的方式提取植被,mDataset为原始影像，
// mNDVIset为NDVI数据，其值在[-1,+1]之间，VegetationSet为提取出的植被数据,th1为低阈值，th2为高阈值
void ExtractVegetation(GDALDataset* mDataset, GDALDataset* mNDVIset, GDALDataset* &VegetationSet, float th1, float th2);

函数文件 Func.cpp

#include "Func.h"void CalcuNDVI(GDALDataset* mDataset, GDALDataset* &mNDVIset)
{// 获取影像的参数int mX = mDataset->GetRasterXSize();int mY = mDataset->GetRasterYSize();int mDataType = mDataset->GetRasterBand(1)->GetRasterDataType();// 将影像中的四个波段取出,分别为NIR，R,G,BGDALRasterBand* mB = mDataset->GetRasterBand(1);GDALRasterBand* mG = mDataset->GetRasterBand(2);GDALRasterBand* mR = mDataset->GetRasterBand(3);GDALRasterBand* mNIR = mDataset->GetRasterBand(4);// 申请缓冲区float* RedBuf = (float*)CPLMalloc(sizeof(float)*mX*mY);float* NIRBuf = (float*)CPLMalloc(sizeof(float)*mX*mY);float* mNDVIband = (float*)CPLMalloc(sizeof(float)*mX*mY);// 从R和NIR波段中取出数据放到缓冲区mR->RasterIO(GF_Read, 0, 0, mX, mY, RedBuf, mX, mY, GDT_UInt32, 0, 0);mNIR->RasterIO(GF_Read, 0, 0, mX, mY, NIRBuf, mX, mY, GDT_UInt32, 0, 0);// 计算NDVIfor (int i = 0; i < mY; i++){for (int j = 0; j < mX; j++){float r = RedBuf[i*mX + j];float nir = NIRBuf[i*mX + j];mNDVIband[i*mX + j] = (nir - r) / (nir + r);}}// 将计算结果写入NDVI的数据集中mNDVIset->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, mX, mY, mNDVIband, mX, mY, GDT_Float32, 0, 0);// 释放内存CPLFree(RedBuf);CPLFree(NIRBuf);CPLFree(mNDVIband);}void ExtractVegetation(GDALDataset * mDataset, GDALDataset * mNDVIset, GDALDataset *& VegetationSet, float th1, float th2)
{// 获取影像参数int mX = mDataset->GetRasterXSize();int mY = mDataset->GetRasterYSize();GDALDataType mImgDataType = mDataset->GetRasterBand(1)->GetRasterDataType();GDALDataType mNDVIDataType = mNDVIset->GetRasterBand(1)->GetRasterDataType();// 获取原始影像的四个波段GDALRasterBand* mBand1 = mDataset->GetRasterBand(1);GDALRasterBand* mBand2 = mDataset->GetRasterBand(2);GDALRasterBand* mBand3 = mDataset->GetRasterBand(3);GDALRasterBand* mBand4 = mDataset->GetRasterBand(4);// 获取NDVI的波段GDALRasterBand* mNDVIBand = mNDVIset->GetRasterBand(1);// 申请缓冲区，用来保存原始波段和NDVI的数据unsigned short* mBandBuf1 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);unsigned short* mBandBuf2 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);unsigned short* mBandBuf3 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);unsigned short* mBandBuf4 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);float* mNDVIBandBuf = (float*)CPLMalloc(mX*mY * mNDVIDataType);// 将各个波段的数据放入对应的缓冲区mBand1->RasterIO(GF_Read, 0, 0, mX, mY, mBandBuf1, mX, mY, mImgDataType, 0, 0);mBand2->RasterIO(GF_Read, 0, 0, mX, mY, mBandBuf2, mX, mY, mImgDataType, 0, 0);mBand3->RasterIO(GF_Read, 0, 0, mX, mY, mBandBuf3, mX, mY, mImgDataType, 0, 0);mBand4->RasterIO(GF_Read, 0, 0, mX, mY, mBandBuf4, mX, mY, mImgDataType, 0, 0);mNDVIBand->RasterIO(GF_Read, 0, 0, mX, mY, mNDVIBandBuf, mX, mY, mNDVIDataType, 0, 0);// 申请缓冲区，用来保存提取出的植被区域unsigned short* mVegBand1 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);unsigned short* mVegBand2 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);unsigned short* mVegBand3 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);unsigned short* mVegBand4 = (unsigned short*)CPLMalloc(mX*mY * mImgDataType);//  通过阈值来提取植被区域for (int i = 0; i < mX*mY; i++){// 若在阈值范围内，赋对应的值，若不在，赋值为0if (mNDVIBandBuf[i] >= th1&&mNDVIBandBuf[i] <= th2){mVegBand1[i] = mBandBuf1[i];mVegBand2[i] = mBandBuf2[i];mVegBand3[i] = mBandBuf3[i];mVegBand4[i] = mBandBuf4[i];}else{mVegBand1[i] = 0;mVegBand2[i] = 0;mVegBand3[i] = 0;mVegBand4[i] = 0;}}// 把提取结果放到数据集中VegetationSet->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, mX, mY, mVegBand1, mX, mY, mImgDataType, 0, 0);VegetationSet->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, mX, mY, mVegBand2, mX, mY, mImgDataType, 0, 0);VegetationSet->GetRasterBand(3)->RasterIO(GF_Write, 0, 0, mX, mY, mVegBand3, mX, mY, mImgDataType, 0, 0);VegetationSet->GetRasterBand(4)->RasterIO(GF_Write, 0, 0, mX, mY, mVegBand4, mX, mY, mImgDataType, 0, 0);// 释放内存CPLFree(mBandBuf1);CPLFree(mBandBuf2);CPLFree(mBandBuf3);CPLFree(mBandBuf4);CPLFree(mNDVIBandBuf);CPLFree(mVegBand1);CPLFree(mVegBand2);CPLFree(mVegBand3);CPLFree(mVegBand4);
}

运行结果

影像数据：

NDVI：

提取植被：

结语

这一个实验中间确实卡了好几次，中间也在弄其他东西，以上的代码里面确实还有一些不合适的地方，但，，，总之结果是跑出来了ㄟ(川.一ㄟ)。后面再慢慢改进吧！主要还是对于C语言的理解不够，很多东西都是边用边学，一知半解的，看官也请多包涵。总之，任重而道远，慢慢学吧~~

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！