基于GDAL计算NDVI

// NDVI.cpp : 定义控制台应用程序的入口点。//#include "stdafx.h"//#include "gdal_priv.h"#include "gdal.h"void CalcNDVI(int rows, int cols, float* nir, float* red, float* ndvi){int n = rows*cols;for (int i = 0; i < n; i ++){if (fabs(nir[i]+red[i]) < 0.1)ndvi[i] = -1;elsendvi[i] = (nir[i]-red[i])/(nir[i]+red[i]);}}int SaveNDVI(int rows, int cols, float* ndvi, const char* ndviFile){GDALDriverH hDriver = GDALGetDriverByName ("HFA");if (hDriver == NULL)return -1;GDALDatasetH hNDVI = GDALCreate(hDriver, ndviFile, cols, rows, 1, GDALDataType::GDT_Float32, NULL);if (hNDVI == NULL)return -1;GDALRasterBandH hNdviBand = GDALGetRasterBand(hNDVI, 1);CPLErr r = GDALRasterIO(hNdviBand, GDALRWFlag::GF_Write, 0, 0,cols, rows, (void*)ndvi, cols, rows, GDALDataType::GDT_Float32,0,0); GDALClose(hNDVI);if (CE_None != r)return -1;return 0;}int _tmain(int argc, _TCHAR* argv[]){GDALAllRegister();//CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");char tmFile[] = "G:\\can_tmr.img"; //数据文件要与头文件放在一起char ndviFile[] = "G:\\ndvi.img";GDALDatasetH hTM = GDALOpen (tmFile, GA_ReadOnly);if (hTM == NULL){printf("error");exit(-1);}int bandCount = GDALGetRasterCount(hTM);GDALRasterBandH hNir = GDALGetRasterBand(hTM, 4);//第四波段为近红外波段;GDALRasterBandH hRed = GDALGetRasterBand(hTM, 3);//第三波段为红波段；int rows = GDALGetRasterYSize(hTM);int cols = GDALGetRasterXSize(hTM);float* nir = new float[rows*cols];float* red = new float[rows*cols];float* ndvi = new float[rows*cols];CPLErr r1 = GDALRasterIO(hNir, GDALRWFlag::GF_Read, 0, 0,cols, rows, (void*)nir, cols, rows, GDALDataType::GDT_Float32,0,0); CPLErr r2 = GDALRasterIO(hRed, GDALRWFlag::GF_Read, 0, 0,cols, rows, (void*)red, cols, rows, GDALDataType::GDT_Float32,0,0); CalcNDVI(rows, cols, nir, red, ndvi);SaveNDVI(rows, cols, ndvi, ndviFile);GDALClose(hTM);GDALDestroyDriverManager();system("pause");return 0;}