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

利用GDAL进行波段合成

admin2个月前 (10-01)遥感94

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

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

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

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

tools.h头文件:

#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实现文件:

#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;
	}
}

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


版权声明:本文由cyhu's essay发布,如需转载请注明出处。

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