使用GDAL对HDF数据进行geoloc校正

在上一篇博客中，大概说了下怎么使用gdal提供的gdalwarp工具来进行校正处理。其实质与envi的glt校正应该是一样的。我把gdalwarp的代码封装了一下，写了一个类来进行geoloc处理。希望对大家有用。

先是头文件，函数的注释很详细，就不多说了。后面的源文件就是从gdalwarp.cpp中摘录出来的，有兴趣的可以看gdalwarp.cpp ，下面的代码只是把这个文件中没有用到的代码删除了，同时有些参数直接写死了。不太清楚的可以直接看原来的代码。

/***************************************************************************
*
* Time: 2013-03-02
* Project: 遥感处理算法 
* Purpose: 地理查找表校正
* Author:  李民录
* Copyright (c) 2013, liminlu0314@gmail.com
* Describe: 地理查找表校正，用于HDF数据的校正算法
*
****************************************************************************/

#ifndef IMAGEGEOLOCWARP_H  
#define IMAGEGEOLOCWARP_H

#include "ImgCore.h"

/**
* \file ImageGeoLocWarp.h 
* @brief 地理查找表校正
*
* 地理查找表校正，主要用于HDF数据的校正。比如海洋卫星数据等，MODIS数据等。
* 一般在该类型的数据中，有两个Latitude和Longitude的32位浮点数的数据集，
* 用来存储数据的经纬度坐标，在校正过程中使用这两个进行处理
*/
class IMGALG_API CImageGeoLocWarp
{
public:
	/**
	* @brief 构造函数，地理查找表校正
	* @param pszSrcFile		原始数据路径
	* @param pszDstFile		结果数据路径
	* @param pszFormat		结果数据格式
	* @param pProcess		进度条指针
	*/
	CImageGeoLocWarp(const char* pszSrcFile, const char* pszDstFile,
		const char* pszFormat = "HFA", CProcessBase* pProcess = NULL);

	/**
	* @brief 析构函数
	*/
	~CImageGeoLocWarp();

	/**
	* @brief 获取输出参数信息
	* @param iHeight			输出图像高度
	* @param iWidth				输出图像宽度
	* @param padfGeoTransform	输出GeoTransform六参数
	* @param padfExtent			输出图像四至范围
	* @return	是否计算成功，以及错误代码
	*/
	int GetSuggestedWarpOutput(int &iHeight, int &iWidth, double *padfGeoTransform = NULL,	double *padfExtent = NULL);

	/**
	* @brief 执行校正算法
	* @param dResX			输出图像横向分辨率，默认为0，表示自动计算输出象元大小
	* @param dResY			输出图像纵向分辨率，默认为0，表示自动计算输出象元大小
	* @param resampling		重采样方式，默认为双线性内插
	* @return	是否计算成功，以及错误代码
	*/
	int DoGeoLocWarp(double dResX = 0.0, double dResY = 0.0, ResampleMethod resampling = RM_Bilinear);

private:
	//输入参数
	const char* m_pszSrcFile;	/*! 原始数据路径 */
	const char* m_pszDstFile;	/*! 结果数据路径 */
	const char* m_pszFormat;	/*! 结果数据格式 */
	CProcessBase* m_pProcess;	/*! 进度条指针 */
	ResampleMethod m_Resampling; /*! 重采样方式 */

	string m_strWkt;	/*! 输出文件投影串，这里只能是WGS84经纬度 */
};

#endif //IMAGEGEOLOCWARP_H

下面是源文件，主要是对gdalwarp.cpp进行了精简，去除了和geoloc无关的代码。

/***************************************************************************
*
* Time: 2013-03-02
* Project: 遥感处理算法
* Purpose: 地理查找表校正
* Author:  李民录
* Copyright (c) 2013, liminlu0314@gmail.com
* Describe: 地理查找表校正，用于HDF数据的校正算法
*
****************************************************************************/

#include "ImageGeoLocWarp.h"//上面的头文件
#include "ImageCommon.h"//这个头文件就是之前的那个使用gdal对数据进行管理的博客（用到的函数是删除临时文件）
#include "gdal_priv.h"
#include "gdalwarper.h"
#include "ogr_srs_api.h"

CImageGeoLocWarp::CImageGeoLocWarp(const char* pszSrcFile, const char* pszDstFile,
					const char* pszFormat, CProcessBase* pProcess)
{
	m_pszSrcFile = pszSrcFile;	/*! 原始数据路径 */
	m_pszDstFile = pszDstFile;	/*! 结果数据路径 */
	m_pszFormat = pszFormat;	/*! 结果数据格式 */
	m_pProcess = pProcess;	/*! 进度条指针 */
	m_Resampling = RM_Bilinear; /*! 重采样方式，默认为双线性 */
	m_strWkt = SRS_WKT_WGS84;	//投影只能是WGS84
}

CImageGeoLocWarp::~CImageGeoLocWarp()
{
}

int CImageGeoLocWarp::GetSuggestedWarpOutput(int &iHeight, int &iWidth, double *padfGeoTransform, double *padfExtent)
{
	if(m_pszSrcFile == NULL || m_pszDstFile == NULL)
		return RE_PARAMERROR;

	GDALAllRegister();

	GDALDatasetH hSrcDS = GDALOpen( m_pszSrcFile, GA_ReadOnly );	//打开文件
	if(hSrcDS == NULL )
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("打开待纠正影像失败！");

		return RE_FILENOTEXIST;
	}

	if(m_strWkt.empty())	//目标投影不存在
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("目标投影不存在！");

		GDALClose( hSrcDS );
		return RE_PARAMERROR;
	}

	// 创建坐标转换参数
	char** papszWarpOptions = NULL;
	papszWarpOptions = CSLSetNameValue( papszWarpOptions, "METHOD", "GEOLOC_ARRAY" );
	papszWarpOptions = CSLSetNameValue( papszWarpOptions, "DST_SRS", m_strWkt.c_str() );

	void *hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszWarpOptions );
	if(hTransformArg == NULL )
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("坐标转换参数错误！");

		CSLDestroy( papszWarpOptions );
		GDALClose( hSrcDS );

		return RE_PARAMERROR;
	}

	GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;

	// 获取目标文件大小以及转换参数信息
	double adfDstGeoTransform[6] = {0};
	double adfExtent[4];
	int nPixels = 0, nLines = 0;

	CPLErr eErr = GDALSuggestedWarpOutput2( hSrcDS, psInfo->pfnTransform, hTransformArg, 
		adfDstGeoTransform, &nPixels, &nLines, adfExtent, 0 );

	iWidth = nPixels;
	iHeight = nLines;
	if(padfGeoTransform != NULL)
		memcpy(padfGeoTransform, adfDstGeoTransform, sizeof(double)*6);
	if(padfExtent != NULL)
		memcpy(padfExtent, adfExtent, sizeof(double)*4);

	GDALDestroyGenImgProjTransformer(hTransformArg);
	CSLDestroy( papszWarpOptions );
	GDALClose( hSrcDS );

	if(eErr != CE_None)
		return RE_PARAMERROR;
	else
		return RE_SUCCESS;
}

int CImageGeoLocWarp::DoGeoLocWarp(double dResX, double dResY, ResampleMethod resampling)
{
	m_Resampling = resampling;

	if(m_pszSrcFile == NULL || m_pszDstFile == NULL)
		return RE_PARAMERROR;

	if(FileIsUsed(m_pszDstFile, m_pProcess))
		return RE_FILEISUESED;

	if(m_pProcess != NULL)
	{
		m_pProcess->ReSetProcess();
		m_pProcess->SetMessage("正在执行影像地理校正...");
	}

	GDALAllRegister();

	GDALDatasetH hSrcDS = GDALOpen( m_pszSrcFile, GA_ReadOnly );	//打开文件
	if(hSrcDS == NULL )
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("打开待纠正影像失败！");

		return RE_FILENOTEXIST;
	}

	if(m_strWkt.empty())	//目标投影不存在
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("目标投影不存在！");

		GDALClose( hSrcDS );
		return RE_PARAMERROR;
	}

	// 创建坐标转换参数
	char** papszWarpOptions = NULL;
	papszWarpOptions = CSLSetNameValue( papszWarpOptions, "METHOD", "GEOLOC_ARRAY" );
	papszWarpOptions = CSLSetNameValue( papszWarpOptions, "DST_SRS", m_strWkt.c_str() );

	void *hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszWarpOptions );
	if(hTransformArg == NULL )
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("坐标转换参数错误！");

		CSLDestroy( papszWarpOptions );
		GDALClose( hSrcDS );

		return RE_PARAMERROR;
	}

	GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;

	// 获取目标文件大小以及转换参数信息
	double adfDstGeoTransform[6] = {0};
	double adfExtent[4];
	int nPixels = 0, nLines = 0;

	CPLErr eErr = GDALSuggestedWarpOutput2( hSrcDS, psInfo->pfnTransform, hTransformArg, 
		adfDstGeoTransform, &nPixels, &nLines, adfExtent, 0 );

	if(eErr != CE_None)
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("获取目标文件大小以及转换参数信息错误！");

		GDALDestroyGenImgProjTransformer(hTransformArg);
		CSLDestroy( papszWarpOptions );
		GDALClose( hSrcDS );

		return RE_PARAMERROR;
	}

	// 如果用户指定了输出图像的分辨率
	if ( dResX != 0.0 || dResY != 0.0 )
	{
		// 如果只指定了一个，使用自动计算的结果
		if ( dResX == 0.0 ) dResX = adfDstGeoTransform[1];
		if ( dResY == 0.0 ) dResY = adfDstGeoTransform[5];

		// 确保分辨率符号正确
		if ( dResX < 0.0 ) dResX = -dResX;
		if ( dResY > 0.0 ) dResY = -dResY;

		// 计算输出图像的范围
		double minX = adfDstGeoTransform[0];
		double maxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
		double maxY = adfDstGeoTransform[3];
		double minY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;

		// 按照用户指定的分辨率来计算图像的输出大小以及范围
		nPixels = ( int )((( maxX - minX ) / dResX ) + 0.5 );
		nLines  = ( int )((( minY - maxY ) / dResY ) + 0.5 );
		adfDstGeoTransform[0] = minX;
		adfDstGeoTransform[3] = maxY;
		adfDstGeoTransform[1] = dResX;
		adfDstGeoTransform[5] = dResY;
	}

	// 获得波段数目
	int nBandCount = GDALGetRasterCount(hSrcDS);
	// 创建输出文件数据类型与源文件相同
	GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));

	// 创建输出结果文件
	GDALDriverH hDriver = GDALGetDriverByName( m_pszFormat);
	if(hDriver == NULL )
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("不支持该种类型数据！");

		GDALDestroyGenImgProjTransformer(hTransformArg);
		CSLDestroy( papszWarpOptions );
		GDALClose( hSrcDS );

		return RE_FILENOTSUPPORT;
	}

	GDALDatasetH  hDstDS = GDALCreate( hDriver, m_pszDstFile, nPixels, nLines, nBandCount, eDT, NULL );
	if(hDstDS == NULL)
	{
		if(m_pProcess != NULL)
			m_pProcess->SetMessage("创建输出文件失败！");

		GDALDestroyGenImgProjTransformer(hTransformArg);
		CSLDestroy( papszWarpOptions );
		GDALClose( hSrcDS );

		return RE_CREATEFAILED;
	}

	// 设置投影等信息
	GDALSetProjection( hDstDS, m_strWkt.c_str() );
	GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

	if (hTransformArg != NULL)
		GDALSetGenImgProjTransformerDstGeoTransform( hTransformArg, adfDstGeoTransform);

	// 逐个波段获取颜色表
	GDALColorTableH hCT;

	for (int i=1; i<=nBandCount; i++)//band从1开始算
	{
		hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS, i) );
		if( hCT != NULL )
			GDALSetRasterColorTable( GDALGetRasterBand(hDstDS, i), hCT );
	}

	// 创建转换选项
	CSLDestroy( papszWarpOptions );	//先释放之前的
	papszWarpOptions = NULL;
	papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");

	GDALTransformerFunc pfnTransformer = NULL;
	pfnTransformer = GDALGenImgProjTransform;	//坐标转换函数

	GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
	psWarpOptions->papszWarpOptions = CSLDuplicate(papszWarpOptions);
	psWarpOptions->hSrcDS = hSrcDS;
	psWarpOptions->hDstDS = hDstDS;
	psWarpOptions->nBandCount = nBandCount;

	// 申请内存
	psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );  
	psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );

	for (int j=0; j<nBandCount; j++)
	{
		psWarpOptions->panSrcBands[j] = j + 1;  
		psWarpOptions->panDstBands[j] = j + 1;
	}

	psWarpOptions->eResampleAlg = (GDALResampleAlg)resampling;
	psWarpOptions->pfnProgress = ALGTermProgress;	//进度条回调函数;  
	psWarpOptions->pProgressArg = m_pProcess;	//进度条指针
	psWarpOptions->papszWarpOptions = papszWarpOptions;

	psWarpOptions->pfnTransformer = pfnTransformer;
	psWarpOptions->pTransformerArg = hTransformArg;

	CPLErr Err = CE_None;
	GDALWarpOperation oOperation;
	if( oOperation.Initialize( psWarpOptions ) == CE_None )
		Err = oOperation.ChunkAndWarpImage( 0, 0, nPixels, nLines);
	else
	{
		CSLDestroy( papszWarpOptions );
		GDALDestroyGenImgProjTransformer( hTransformArg );

		GDALClose( hDstDS );
		GDALClose( hSrcDS );

		RasterDelete(m_pszDstFile);	//删除结果数据
		return RE_FAILED;
	}

	CSLDestroy( papszWarpOptions );
	GDALDestroyGenImgProjTransformer( hTransformArg );

	GDALClose( hDstDS );
	GDALClose( hSrcDS );

	if (m_pProcess && !m_pProcess->m_bIsContinue)
	{
		RasterDelete(m_pszDstFile);	//删除结果数据
		return RE_USERCANCEL;
	}

	if (Err != CE_None)
		return RE_FAILED;

	return RE_SUCCESS;
}

调用代码如下：

#include <stdio.h>
#include <stdlib.h>

#include <boost/progress.hpp>	//boost计时函数

int main()
{
	int iRev = RE_SUCCESS;

	CConsoleProcess *pPro = new CConsoleProcess();
	progress_timer *pTime = new progress_timer();  // 开始计时

	// 对HDF数据中的一个字数据进行处理，下面分别是子数据集的路径和输出文件路径
	const char* pszSrc = "HDF4_SDS:UNKNOWN:\"D:\\Data\\hdf\\H1BCLR101106020418636.L2C.HDF\":12";
	const char* pszDst = "D:\\Data\\hdf\\H1BCLR101106020418636_12.tif";

	CImageGeoLocWarp warp(pszSrc, pszDst, "HFA", pPro);
	iRev = warp.DoGeoLocWarp(0.0, 0.0, RM_NearestNeighbour);

	if(iRev == RE_SUCCESS)
		printf("success!\n");
	else
		printf("failed!\n");

	delete pPro;
	delete pTime;
	::system("pause");
	return 0;
}

最后再说一下，代码中的返回值，RE_开头的去我之前的博客中去找，就是一些常量，用于标记返回的信息。还有那个进度条的类，也在我之前的博客都有说明的。上面的代码中用到的第三方库就是BOOST和GDAL，算法中只有GDAL一个，BOOST库只在main函数中用来计时用的。其他的在我之前的博客都有说明的。

不要指望把我的代码拿过去直接编译就能过去，那样的话，我觉得对你也没有啥进步。如果一个人看懂我的代码的话，那些不认识的函数也就大概能猜出是个啥意思了，自己写一个也很简单的。学习别人的代码，不要指望直接就能用，关键是还要看别人的思路，然后结合自己的项目（工程）来进行改造来达到自己的目的，这才是一种对自己有利的学习方式。废话有点多，觉得说的没道理的后面这些话就当没有吧。