当前位置：首页 > 遥感 > 正文内容

利用GDAL进行波段合成

admin3年前 (2021-10-01)遥感11047

顺手写的一个小工具，给出头文件和实现

主要用了GDALDataset类下面的RasterIO函数，没有用GDALBand下面的RasterIO

另外实现了一个定时器类，实际上就是之前写的一个<http://lovedm.club/?id=98>

输出文件格式为固定的GTiff，函数的第三个参数是GTiff创建选项，相关选项参考<https://gdal.org/drivers/raster/gtiff.html>

tools.h头文件：

``` c++
#pragma once
#include <gdal_priv.h>
#include <vector>
#include <chrono>
#include <string>

using namespace std;
using namespace std::chrono;

static class tools
{
public:
	/// <summary>
	/// 计时器类
	/// </summary>
	class HightPrecisionTimeStamp
	{
	public:
		// 刷新当前时间点
		void Update();
		// 获取时间秒
		double GetElapsedTimeInSecond();
		// 获取时间毫秒
		double GetElapsedTimeInMilliSec();
		// 获取时间微秒
		long long GetElapsedTimeInMicroSec();
	protected:
		time_point<high_resolution_clock> m_begin; //时间开始
	};
	/// <summary>
	/// 波段合成
	/// 输出影像类型为GTIFF
	/// </summary>
	/// <param name="paths">单个波段的全路径</param>
	/// <param name="outPath">输出全路径</param>
	/// <param name="papszOptions">tiff创建选项</param>
	/// <returns>是否成功</returns>
	static bool LayerStack(vector<string> paths, string outPath, char** papszOptions = NULL);

};
```

tools.cpp实现文件：

``` c++
#include "tools.h"

void tools::HightPrecisionTimeStamp::Update()
{
	m_begin = high_resolution_clock::now();
}
double tools::HightPrecisionTimeStamp::GetElapsedTimeInSecond()
{
	return GetElapsedTimeInMicroSec() * 0.000001;
}
double tools::HightPrecisionTimeStamp::GetElapsedTimeInMilliSec()
{
	return GetElapsedTimeInMicroSec() * 0.001;
}
long long tools::HightPrecisionTimeStamp::GetElapsedTimeInMicroSec()
{
	return duration_cast<microseconds>(high_resolution_clock::now() - m_begin).count();
}

bool tools::LayerStack(vector<string> paths, string outPath, char** papszOptions)
{
	GDALAllRegister();
	vector<GDALDataset*>dataSets;
	try
	{

for (int i = 0; i < paths.size(); i++)
		{
			GDALDataset* tempDataset = (GDALDataset*)GDALOpen(paths[i].c_str(), GA_ReadOnly);
			dataSets.push_back(tempDataset);
		}

int rasterXSize = dataSets[0]->GetRasterXSize();
		int rasterYSize = dataSets[0]->GetRasterYSize();
		double geotransform[6];
		dataSets[0]->GetGeoTransform(geotransform);
		const char* project = dataSets[0]->GetProjectionRef();
		GDALDataType dataType = dataSets[0]->GetRasterBand(1)->GetRasterDataType();

if (papszOptions == NULL)
		{
			papszOptions = CSLSetNameValue(papszOptions, "TILED", "YES");
			papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "LZW");
			papszOptions = CSLSetNameValue(papszOptions, "NUM_THREADS", "ALL_CPUS");
			papszOptions = CSLSetNameValue(papszOptions, "INTERLEAVE", "PIXEL");//BAND是BSQ格式，PIXEL是BIP格式，默认是PIXEL
		}

GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTIFF");
		GDALDataset* outDataset = poDriver->Create(outPath.c_str(), rasterXSize, rasterYSize, paths.size(), dataType, papszOptions);

outDataset->SetGeoTransform(geotransform);
		outDataset->SetProjection(project);

void* pdata = NULL;
		if (dataType == GDT_Byte)
		{
			pdata = CPLMalloc(sizeof(unsigned char) * rasterXSize * rasterYSize);
		}
		else if (dataType == GDT_Int16)
		{
			pdata = CPLMalloc(sizeof(int16_t) * rasterXSize * rasterYSize);
		}
		else if (dataType == GDT_UInt16)
		{
			pdata = CPLMalloc(sizeof(uint16_t) * rasterXSize * rasterYSize);
		}
		else if (dataType == GDT_Int32)
		{
			pdata = CPLMalloc(sizeof(int32_t) * rasterXSize * rasterYSize);
		}
		else if (dataType == GDT_UInt32)
		{
			pdata = CPLMalloc(sizeof(uint32_t) * rasterXSize * rasterYSize);
		}
		else if (dataType == GDT_Float32)
		{
			pdata = CPLMalloc(sizeof(float) * rasterXSize * rasterYSize);
		}
		else if (dataType == GDT_Float64)
		{
			pdata = CPLMalloc(sizeof(double) * rasterXSize * rasterYSize);
		}
		//unsigned short* pdata = new unsigned short[xSize * ySize];
		//unsigned short* pdata = (unsigned short*)CPLMalloc(sizeof(unsigned short) * rasterXSize * rasterYSize);

int bandlist1[1] = { 1 };
		int bandlist2[1] = { 1 };
		for (int i = 0; i < paths.size(); i++)
		{
			bandlist2[0] = i + 1;
			CPLErr err1 = dataSets[i]->RasterIO(GF_Read, 0, 0, rasterXSize, rasterYSize, pdata, rasterXSize, rasterYSize, dataType, 1, bandlist1, 0, 0, 0, NULL);
			CPLErr err2 = outDataset->RasterIO(GF_Write, 0, 0, rasterXSize, rasterYSize, pdata, rasterXSize, rasterYSize, dataType, 1, bandlist2, 0, 0, 0, NULL);
			GDALClose((GDALDatasetH)dataSets[i]);
			if (err1 == CE_Failure || err2 == CE_Failure)
			{
				CPLFree(pdata);
				GDALClose((GDALDatasetH)outDataset);
				return false;
			}
		}

CPLFree(pdata);
		GDALClose((GDALDatasetH)outDataset);
		return true;
	}
	catch (exception e)
	{
		return false;
	}
}
```

**随手写的，有问题请发邮件告知。**

扫描二维码推送至手机访问。

本文链接：http://lovedm.club/?id=137

标签: C++

分享给朋友：

返回列表

上一篇：GDAL-3.0.1+PROJ-6.2.1+SQLite+HDF4+HDF5+VS2019编译

下一篇：Python embed版本安装使用

“利用GDAL进行波段合成” 的相关文章

个人学习笔记

利用GDAL进行波段合成

“利用GDAL进行波段合成” 的相关文章

太阳天顶角、太阳高度角、太阳方位角

Python embed版本安装使用

国产卫星绝对辐射定标系数（2008——2022）

Python GDAL In GDAL 4.0, exceptions will be enabled by default

鲁ICP备18000987号-1 鲁公网安备 37148202000241号

Powered By Z-BlogPHP. Theme by TOYEAN.