注册 登录
  • 注册时,本站名称为:RGB空格3S博客,注意中间的空格。
  • 2018.10.14本人利用业余时间接遥感数据处理与编程小活,欢迎骚扰,QQ:853060844
  • 2018.8.14本站全面接入google广告
  • 本人闲暇时间(周末及其他闲暇时间)接遥感图像处理与IDL编程小活,欢迎骚扰,qq:853060844
  • 2017.2.14今天收到45条恶意评论(全是外文),故评论时请填写必要信息,匿名评论全部拉黑,迫不得已而为之
  • 2017.1.27,2017年春节及至,我谨代表本人祝大家新春快乐,本人年终总结文章请访问:2016年终总结
  • 为防止恶意转载,本站全面禁止复制,并添加图片水印:RGB 3S博客www.ixxin.cn。
  • 本站正式更名为RGB 3S博客,本站将撤消所有非3S内容,其将转移到新博客江湖时代

Landsat8 GDAL C++读取接口

C/C++ admin 488次浏览 已收录 2个评论
[隐藏]

前言

Landsat8数据是最常用的免费数据,其具有9个可见光短波红外波段,2个热红外波段,可以反演多个地表、大气参数,包括大气气溶胶、地表温度。然而其下载文件解压后,由多个tif组成多个文件,而不是一个tif中有多个波段,这对GDAL读取造成了一定的困扰,但是莫慌,这里写了一个GDAL接口用来读取Landsat8数据:

具体实现

首先是一个头文件“Landsat8Dataset.h”

#pragma once
#include "C:\warmerda\x64\include\gdal_priv.h"
#include <vector>
#include <string>
#include "gdal_utils.h"
#include "cplkeywordparser.h"
#include <iostream>
#include "gdal_pam.h"
using namespace std;

enum L8DTS
{
	Pan,
	Mul,
	QA,
	Them
};
class Landsat8Dataset
{
public:
	Landsat8Dataset();
	~Landsat8Dataset();
	int nBands = 0;
	GDALDataset* poVRTDS;
	GDALDatasetH* GetDatasetList();
	char** FileList = NULL;
	GDALDataset* Open(const char* pMTLFile,L8DTS DT);
	void Close();
protected:
	GDALDatasetH* DatasetList = NULL;
};

具体实现文件“Landsat8Dataset.cpp”

#include "Landsat8Dataset.h"

Landsat8Dataset::Landsat8Dataset()
{
	poVRTDS = NULL;
}


Landsat8Dataset::~Landsat8Dataset()
{
	Close();
}

GDALDataset* Landsat8Dataset::Open(const char* pMTLFile,L8DTS DT = L8DTS::Mul)
{
	/*初始化参数*/
	double dfCellSize = 0.0;
	int nRasterXSize = 0, nRasterYSize = 0;
	char **papszRasterName = NULL;
	int nSubDataset = 1;
	switch (DT)
	{
	case Pan:
		nSubDataset = 0;
		break;
	case Mul:
		nSubDataset = 1;
		break;
	case QA:
		nSubDataset = 4;
		break;
	case Them:
		nSubDataset = 2;
		break;
	default:
		break;
	}
	
	GDALAllRegister();
	/*读取MTL文件*/
	VSILFILE *fp = VSIFOpenL(pMTLFile, "r");
	if (fp == NULL)
	{
		return NULL;
	}

	CPLKeywordParser oParser;
	if (!oParser.Ingest(fp))
	{
		VSIFCloseL(fp);
		return NULL;
	}

	VSIFCloseL(fp);

	char** papszMetaInfo = oParser.GetAllKeywords();

	if (nSubDataset == 0)
	{
		dfCellSize = CPLAtof(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PROJECTION_PARAMETERS.GRID_CELL_SIZE_PANCHROMATIC", "0"));
		nRasterXSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.PANCHROMATIC_SAMPLES", "0"));
		nRasterYSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.PANCHROMATIC_LINES", "0"));
		papszRasterName = CSLAddString(papszRasterName, "BAND_8");
	}
	else if (nSubDataset == 1)
	{
		dfCellSize = CPLAtof(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PROJECTION_PARAMETERS.GRID_CELL_SIZE_REFLECTIVE", "0"));
		nRasterXSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.REFLECTIVE_SAMPLES", "0"));
		nRasterYSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.REFLECTIVE_LINES", "0"));
		papszRasterName = CSLAddString(papszRasterName, "BAND_1");
		papszRasterName = CSLAddString(papszRasterName, "BAND_2");
		papszRasterName = CSLAddString(papszRasterName, "BAND_3");
		papszRasterName = CSLAddString(papszRasterName, "BAND_4");
		papszRasterName = CSLAddString(papszRasterName, "BAND_5");
		papszRasterName = CSLAddString(papszRasterName, "BAND_6");
		papszRasterName = CSLAddString(papszRasterName, "BAND_7");
		papszRasterName = CSLAddString(papszRasterName, "BAND_9");
	}
	else if (nSubDataset == 2)
	{
		dfCellSize = CPLAtof(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PROJECTION_PARAMETERS.GRID_CELL_SIZE_THERMAL", "0"));
		nRasterXSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.THERMAL_SAMPLES", "0"));
		nRasterYSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.THERMAL_LINES", "0"));
		papszRasterName = CSLAddString(papszRasterName, "BAND_10");
		papszRasterName = CSLAddString(papszRasterName, "BAND_11");
	}
	else
	{
		dfCellSize = CPLAtof(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PROJECTION_PARAMETERS.GRID_CELL_SIZE_THERMAL", "0"));
		nRasterXSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.THERMAL_SAMPLES", "0"));
		nRasterYSize = atoi(CSLFetchNameValueDef(papszMetaInfo, "L1_METADATA_FILE.PRODUCT_METADATA.THERMAL_LINES", "0"));
		papszRasterName = CSLAddString(papszRasterName, "BAND_QUALITY");
	}

	if (!GDALCheckDatasetDimensions(nRasterXSize, nRasterYSize))
	{
		CSLDestroy(papszRasterName);
		return NULL;
	}

	int nBandCount = CSLCount(papszRasterName);
	nBands = nBandCount;
	//GDALDatasetH* DatasetS = new GDALDatasetH[nBandCount]();
	DatasetList = new GDALDatasetH[nBandCount]();
	std::string sTmpStr;
	CPLString ostrFileName;
	const char *pszFilename;
	const char* pszSubFile;
	const char* osDirname = CPLGetDirname(pMTLFile);
	for (int iBand = 0; iBand < nBandCount; iBand++)
	{
		//ostrFileName.Printf("L1_METADATA_FILE.PRODUCT_METADATA.FILE_NAME_%s", papszRasterName[iBand]);
		sTmpStr = papszRasterName[iBand];
		ostrFileName = "L1_METADATA_FILE.PRODUCT_METADATA.FILE_NAME_" + sTmpStr;
		pszFilename = CSLFetchNameValue(papszMetaInfo, ostrFileName.c_str());
		if (pszFilename == NULL)
		{
			CPLError(CE_Failure, CPLE_AppDefined, "Missing %s in .TXT file.", ostrFileName.c_str());
			CSLDestroy(papszRasterName);
			return NULL;
		}

		// trim double quotes. 
		if (pszFilename[0] == '"')
			pszFilename++;
		if (pszFilename[strlen(pszFilename) - 1] == '"')
			((char *)pszFilename)[strlen(pszFilename) - 1] = '\0';

		pszSubFile = CPLFormFilename(osDirname, pszFilename, NULL);
		FileList = CSLAddString(FileList, pszSubFile);
		GDALDatasetH poBandDS = GDALOpen(pszSubFile, GA_ReadOnly);
		//DatasetS[iBand] = poBandDS;
		DatasetList[iBand] = poBandDS;
	}
	char** argv = NULL;
	argv = CSLAddString(argv, "-separate");
	argv = CSLAddString(argv, "-srcnodata");
	argv = CSLAddString(argv, "0");
	argv = CSLAddString(argv, "-vrtnodata");
	argv = CSLAddString(argv, "0");
	GDALBuildVRTOptions* VRTOptions = GDALBuildVRTOptionsNew(argv, NULL);

	GDALDataset* VRTDt = (GDALDataset*)GDALBuildVRT(NULL, nBandCount, DatasetList, NULL, VRTOptions, NULL);
	if (VRTDt == NULL)
	{
		std::cout << "VRT Build Error" << endl;
		GDALBuildVRTOptionsFree(VRTOptions);
		CSLDestroy(papszRasterName);
		for (int i = 0; i < nBandCount; i++)
		{
			GDALClose(DatasetList[i]);
		}
		GDALClose((GDALDatasetH)VRTDt);
		return NULL;

	}
	GDALBuildVRTOptionsFree(VRTOptions);
	CSLDestroy(papszRasterName);
	VRTDt->SetMetadata(papszMetaInfo);
	poVRTDS = VRTDt;
	return VRTDt;
}
void Landsat8Dataset::Close()
{
	if (this->poVRTDS != NULL)
	{
		GDALClose(GDALDatasetH(this->poVRTDS));
	}
	if (this->DatasetList != NULL)
	{
		for (int i = 0; i < nBands; i++)
		{
			if (this->DatasetList[i] == NULL)
				continue;
			else
				GDALClose(this->DatasetList[i]);
		}
	}
	if (this->FileList != NULL)
	{
		CSLDestroy(this->FileList);
	}
	/*归为NULL*/
	this->poVRTDS = NULL;
	this->DatasetList = NULL;
	this->FileList = NULL;
}

读取示例:

#include <iostream>
#include "Landsat8Dataset.h"

using namespace std;
int ReadLandsat8()
{

	const char* pMTLFile = "E:\\Landsat8\\LC81240332013107LGN01\\LC81240332013107LGN01_MTL.txt";
	const char* pOLIFile = "E:\\Landsat8\\OLI.tif";
	const char* pPANFile = "E:\\Landsat8\\PAN.tif";

	Landsat8Dataset* Landsat8 = new Landsat8Dataset();
	GDALDataset* GDALLandsatDt = Landsat8->Open(pMTLFile,L8DTS::Pan);

	int nbands = Landsat8->nBands;

	char** FileList = Landsat8->FileList;

	int nFile = CSLCount(FileList);

	for (int f = 0; f < nFile; f++)
	{
		cout << FileList[f] << endl;
	}
	GDALDriver* pDstDrive = (GDALDriver*)GDALGetDriverByName("ENVI");
	GDALDataset* pDstDs = pDstDrive->CreateCopy(pPANFile, GDALLandsatDt, FALSE, NULL,(GDALProgressFunc)GDALTermProgress, NULL);
	GDALClose(GDALDatasetH(pDstDs));
	/*关闭*/
	Landsat8->Close();

	if (Landsat8->poVRTDS != NULL)
	{
		cout << 1 << endl;
	}

	delete Landsat8;

	cout << "Success" << endl;
	return 0;
}

结尾:
通过此接口就可以直接读取Landsat8 MTL头文件了,是不是很方便?具体实现可以看一下第三段示例代码哦。源代码下载链接:https://download.csdn.net/download/wudixinxin/10636979
,支持博主的点击一下网站头部与底部的广告哦,扫一下支付宝红包。

如有问题,请发邮件到:s_xxin@qq.com


xxin blog , 版权所有丨如未注明 , 均为原创丨本网站采用BY-NC-SA协议进行授权 , 转载请注明Landsat8 GDAL C++读取接口
喜欢 (0)
支付宝[853060844@qq.com]
分享 (0)
admin
关于作者:
坐标山科大遥感系小鲜肉一枚。
发表我的评论
取消评论
表情 贴图 加粗 删除线 居中 斜体 签到

Hi,您需要填写昵称和邮箱!

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址
(2)个小伙伴在吐槽
  1. 谢谢分享,旁边的仓鼠很好玩,喂了5分中,哈哈哈
    Page2018-09-11 16:15 回复 Windows 7 | Chrome 63.0.3239.132
    • admin
      不要点它了,帮我点点广告 :grin:
      admin2018-09-11 17:04 回复 Windows 10 | Chrome 67.0.3396.99