基于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;
		else
			ndvi[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;
}