使用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函数中用来计时用的。其他的在我之前的博客都有说明的。不要指望把我的代码拿过去直接编译就能过去，那样的话，我觉得对你也没有啥进步。如果一个人看懂我的代码的话，那些不认识的函数也就大概能猜出是个啥意思了，自己写一个也很简单的。学习别人的代码，不要指望直接就能用，关键是还要看别人的思路，然后结合自己的项目（工程）来进行改造来达到自己的目的，这才是一种对自己有利的学习方式。废话有点多，觉得说的没道理的后面这些话就当没有吧。
浙公网安备 33010602011771号