python 红外光谱曲线红外光谱数据处理技巧

转载

mob64ca140c3859 2024-07-03 22:16:09

文章标签 python 红外光谱曲线 c++ 图像处理 gdal #include 文章分类 Python 后端开发

文章目录

一、实验任务与内容
二、实验原理
三、编程思路
四、实验成果图
五、源代码

ImageProcessing.h
ImageProcessing.cpp
Main.cpp

一、实验任务与内容

利用VC++编写计算植被指数的程序

（1）采用参数传递的方式，并使用一个变量输出结果

（2）要有说明部分，说明输入参数的格式，说明输出数据的方式

（3）将计算结果用文件保存下来，命名为vegetable.tif

二、实验原理

众所周知，近红外（NIR）的波段为0.78到3微米，而红光（R）的波段为0.625到0.74微米。如图1，在植物的反射波谱特性曲线中，近红外和红光波段中间有较大的落差，这是植物在反射波谱特性曲线中特有的性质，利用这个性质可以在遥感图像中提取植被。

python 红外光谱曲线红外光谱数据处理技巧_#include

图1. 植物的反射波谱特性曲线

定义归一化植被指数（NDVI）为：
$python 红外光谱曲线红外光谱数据处理技巧_图像处理_02$
NDVI越大，表示近红外和红外波段之间的落差越大，则越有可能是植被。通过计算遥感影像的归一化植被指数，再均匀地映射到[0,255]，遥感影像上的像元越白，则越有可能是植被。

三、编程思路

首先创建ImageProcessing类，有六个私有成员，分别为：Dataset，存储读入影像的数据集；xSize、ySize、imageSize、bandNum、dataType分别存储读入影像的列数、行数、像元数、波段数、DN值的数据类型。

private:
	GDALDataset* Dataset;
	int xSize;
	int ySize;
	int imageSize;
	int bandNum;
	GDALDataType dataType;

然后添加公有成员，首先是构造函数和析构函数。构造函数通过传入指向待处理图像路径的字符指针来构造实例。

ImageProcessing(const char*);
~ImageProcessing();

然后创建两个函数calcuNDVI和saveImage分别用来计算NDVI值和保存影像。calcuNDVI函数传入近红外的波段数和红外的波段数，返回一个一维数组，其中的元素是将所有像元的NDVI值由[-1,1]平均映射到[0,255]后的值。saveImage函数传入保存的路径和calcuNDVI返回的一维数组。

unsigned char* calcuNDVI(int, int) const;
void saveImage(const char*, unsigned char*);

在计算NDVI之前，先创建一个函数NDVI2DN，这个函数可以将[-1,1]的值均匀地映射到[0,255]

int NDVI2DN(float);

首先创建四个一维数组，分别用来存储近红外波段的DN值、红外波段的DN值、[-1,1]的NDVI、[0,255]的NDVI

unsigned char* NIRBuffer = new unsigned char[imageSize * sizeof(dataType)];
unsigned char* RBuffer = new unsigned char[imageSize * sizeof(dataType)];
float* NDVIFloatBuffer = new float[imageSize];
unsigned char* NDVIBuffer = new unsigned char[imageSize * sizeof(dataType)];

然后用RasterIO函数将红外波段和近红外波段读入数组

GDALRasterBand* NIRBand = Dataset->GetRasterBand(NIRBandnum);
NIRBand->RasterIO(GF_Read, 0, 0, xSize, ySize, NIRBuffer, xSize, ySize, dataType, 0, 0);
GDALRasterBand* RBand = Dataset->GetRasterBand(RBandnum);
RBand->RasterIO(GF_Read, 0, 0, xSize, ySize, RBuffer, xSize, ySize, dataType, 0, 0);

然后逐像元计算NDVI，然后映射到[0,255]，最后存入NDVIBuffer，删除堆内存中的数组，将NDVIBuffer返回，即完成了遥感影像植被指数的计算。

for(int i = 0; i < ySize; i++) {
	for(int j = 0; j < xSize; j++) {
		int index = xSize * i + j;
		if(int(NIRBuffer[index] + RBuffer[index]) == 0) {
			NDVIFloatBuffer[index] = 0;
		}
		else {
			NDVIFloatBuffer[index] = float(NIRBuffer[index] - RBuffer[index]) / float(NIRBuffer[index] + RBuffer[index]);
		}
		NDVIBuffer[index] = NDVI2DN(NDVIFloatBuffer[index]);
	}
}

delete[] NDVIFloatBuffer;
delete[] NIRBuffer;
delete[] RBuffer;
return NDVIBuffer;

最后将影像保存即可。

使用示例：

首先注册GDAL，然后定义待处理图片路径、保存路径。接着创建ImageProcessing类的实例，调用saveImage方法，在saveImage方法中传出参数，最后删除实例释放内存。

int main() {
	GDALAllRegister();

	const char* imagepath = "..\\image\\test2.tif";
	const char* savepath = "..\\image\\vegetable.tif";

	ImageProcessing* imageProcessing = new ImageProcessing(imagepath);
	//波段4是近红外，波段3是红外
	imageProcessing->saveImage(savepath, imageProcessing->calcuNDVI(4,3));
	delete imageProcessing;
	return 0;
}

四、实验成果图

归一化植被指数均匀映射到[0,255]之后输出的遥感影像如图2所示：

python 红外光谱曲线红外光谱数据处理技巧_c++_03

图2. 实验成果

图中红色实线圈出来的地方都极有可能有植被存在。

五、源代码

ImageProcessing.h

#pragma once
#include<iostream>
#include<gdal.h>
#include<gdal_priv.h>

using namespace std;

class ImageProcessing {
private:
	GDALDataset* Dataset;
	int xSize;
	int ySize;
	int imageSize;
	int bandNum;
	GDALDataType dataType;

public:
	ImageProcessing(const char*);
	~ImageProcessing();
	int getBandNum() const;
	int getImageSize() const;
	void calcuNDVI(const char*, int, int) const;
};

ImageProcessing.cpp

#include<iostream>
#include"ImageProcessing.h"
#include<gdal.h>
#include<gdal_priv.h>
#include<assert.h>

using namespace std;

int NDVI2DN(float);

ImageProcessing::ImageProcessing(const char* path) {
	this->Dataset = (GDALDataset*)GDALOpen(path, GA_ReadOnly);
	this->xSize = Dataset->GetRasterXSize();
	this->ySize = Dataset->GetRasterYSize();
	this->imageSize = this->xSize * this->ySize;
	this->bandNum = Dataset->GetRasterCount();
	this->dataType = Dataset->GetRasterBand(1)->GetRasterDataType();
}

ImageProcessing::~ImageProcessing() {}

int ImageProcessing::getBandNum() const {
	return this->bandNum;
}

int ImageProcessing::getImageSize() const {
	return this->imageSize;
}

void ImageProcessing::calcuNDVI(const char* savepath, int NIRBandnum, int RBandnum) const {
	unsigned char* NIRBuffer = new unsigned char[imageSize * sizeof(dataType)];
	unsigned char* RBuffer = new unsigned char[imageSize * sizeof(dataType)];

	float* NDVIFloatBuffer = new float[imageSize];
	unsigned char* NDVIBuffer = new unsigned char[imageSize * sizeof(dataType)];

	GDALRasterBand* NIRBand = Dataset->GetRasterBand(NIRBandnum);
	NIRBand->RasterIO(GF_Read, 0, 0, xSize, ySize, NIRBuffer, xSize, ySize, dataType, 0, 0);
	GDALRasterBand* RBand = Dataset->GetRasterBand(RBandnum);
	RBand->RasterIO(GF_Read, 0, 0, xSize, ySize, RBuffer, xSize, ySize, dataType, 0, 0);

	for(int i = 0; i < ySize; i++) {
		for(int j = 0; j < xSize; j++) {
			int index = xSize * i + j;
			if(int(NIRBuffer[index] + RBuffer[index]) == 0) {
				NDVIFloatBuffer[index] = 0;
			}
			else {
				NDVIFloatBuffer[index] = float(NIRBuffer[index] - RBuffer[index]) / float(NIRBuffer[index] + RBuffer[index]);
			}
			NDVIBuffer[index] = NDVI2DN(NDVIFloatBuffer[index]);
		}
	}

	GDALDriver* Driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	int BandMap[1] = { 1 };
	char** papszOption = nullptr;
	papszOption = CSLSetNameValue(papszOption, "INTERLEAVE", "BAND");
	GDALDataset* saveDataset = Driver->Create(savepath, this->xSize, this->ySize, 1, this->dataType, papszOption);

	if(!saveDataset) {
		assert(!saveDataset);
	}

	saveDataset->RasterIO(GF_Write, 0, 0, this->xSize, this->ySize, NDVIBuffer, this->xSize, this->ySize, this->dataType, 1, BandMap, 0, 0, 0);

	GDALClose(Dataset);
	GDALClose(saveDataset);

	delete[] NDVIFloatBuffer;
	delete[] NDVIBuffer;
	delete[] NIRBuffer;
	delete[] RBuffer;
}

int NDVI2DN(float NDVI) {
	NDVI += 1;
	int DN;
	if(NDVI == 1) {
		DN = 127;
	}
	else {
		NDVI = NDVI / 2 * 255.0;
		DN = int(NDVI);
	}
	return DN;
}

Main.cpp

#include<iostream>
#include<gdal.h>
#include<gdal_priv.h>
#include"ImageProcessing.h"
using namespace std;

int main() {
	GDALAllRegister();

	const char* imagepath = "..\\image\\test2.tif";
	const char* savepath = "..\\image\\vegetable.tif";

	ImageProcessing* imageProcessing = new ImageProcessing(imagepath);
	//cout << imageProcessing->getBandNum() << endl;
	//cout << imageProcessing->getImageSize();

	//波段4是近红外，波段3是红外
	imageProcessing->calcuNDVI(savepath, 4, 3);
	delete imageProcessing;
	return 0;
}

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。