前言
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++读取接口！
