C++中利用GDAL读取和写入遥感图像，以及和cv::Mat的相互转化

GDAL是一个功能强大的库，但是本身的数据处理算法颇少。常常需要结合OpenCV进行处理。这中间需要进行读取、转化、处理、写出等过程。今天就出于助人为乐的心，向大家展示这些代码。

我的C++水平完全是初学级别的，代码只能供大家参考，一乐。

主函数：

int main()
{	
	//-----------------test_save----------------
	const char* filepath = "F:\\codes\\DIP_C_code\\1\\before.img";
	double geotransform[6]; // 定义地理参考
	char prj[1000] = {""};//用于存储投影
	std::vector<cv::Mat> images = openrsimagewithparameters(filepath, geotransform, prj);
	//输出一下投影信息观看一下
	std::cout << prj << std::endl;
	//好，这之后可以进行某些处理，然后开始进行保存
	const char* savepath = "F:\\codes\\DIP_C_code\\task\\out\\test.tif";
	rsimagesave(images, savepath, prj, geotransform);
	//------------------------------------------------------------
	images.clear();
	return 0;
}

其中用于打开遥感影像并返回std::vectorcv::Mat的函数：该函数可以把投影参数和地理参考的参数都记录下来，用于后面的保存。

std::vector<cv::Mat> openrsimagewithparameters(const char* filepath, double *geotransform,char* prj) {
	GDALDataset *poDataset;
	GDALAllRegister();
	poDataset = (GDALDataset*)GDALOpen(filepath, GA_ReadOnly);
	if (poDataset == NULL) {
		std::cout << "Error:Open Image Failed." << std::endl;
	}
	//some parameters about the dataset.
	int width, height, num_bands;
	width = poDataset->GetRasterXSize();
	height = poDataset->GetRasterYSize();
	num_bands = poDataset->GetRasterCount();
	poDataset->GetGeoTransform(geotransform);
	const char* const_prj = poDataset->GetProjectionRef();
	//std::cout << poDataset->GetProjectionRef() << std::endl;
	char* str= nullptr;
	str = const_cast<char*>(const_prj);
	int length = std::strlen(const_prj);
	for (int i = 0; i < length; i++) {
		prj[i] = str[i];
	}
	//std::cout << num_bands << std::endl;
	std::vector<cv::Mat> images;

	for (int i = 1; i < num_bands + 1; i++) {
		//read, band start from 1.
		GDALRasterBand *poBand;
		float *data_of_band = new float[width * height];
		poBand = poDataset->GetRasterBand(i);
		GDALDataType datatype = poBand->GetRasterDataType();
		poBand->RasterIO(GF_Read, 0, 0, width, height, data_of_band,
			width, height, datatype, 0, 0);//读到缓存
		//new a Mat
		cv::Mat band_img(height, width, CV_8UC1, data_of_band);
		images.push_back(band_img.clone());
	}
	GDALClose(poDataset);
	return images;
}

下面展示rsimagesave函数的代码：

void rsimagesave(std::vector<cv::Mat> images, const char* savefilepath, char* prj, double* geotransform) {
	int channels = images.size();
	int height = images[0].rows;
	int width = images[0].cols;
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset *poDataset;
	poDataset = poDriver->Create(savefilepath, width, height, channels, GDT_Byte, NULL);
	const char* projection = prj;
	poDataset->SetGeoTransform(geotransform);
	poDataset->SetProjection(projection);

	for (int channel = 0; channel < channels; channel++) {
		 GDALRasterBand *poBand = poDataset->GetRasterBand(channel+1);
		 uchar *banddata = new uchar[height * width];
		 cv::Mat strechimage;
		 strechimage = images[channel].clone().reshape(1, 1);
		 for (int col = 0; col < height*width; col++) {
			 banddata[col] = strechimage.at<uchar>(0, col);
		 }
		 poBand->RasterIO(GF_Write, 0, 0, width, height, banddata, width, height, GDT_Byte, 0, 0);
		 std::cout << "-->";
	}
	std::cout << "Finish!" << std::endl;
	GDALClose(poDataset);
}