不知道该怎么详细地写,直接上代码算了! (开发环境的搭建参考我的博文: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; }
本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。
我来说两句