关于GDAL打开hfa大文件的问题[转]

今天在使用GDAL打开大的img文件的时候，（这里所谓的大文件指的是img文件太大，会将数据文件存放到ige格式raw文件中）。在讲img文件和ige文件重命名后，使用GDAL打开文件后，只能读取到文件信息，不能读取图像的数据文件。仔细查看GDAL源代码发现，在img文件中记录了对应的ige文件的名称，重命名后img文件中的记录ige文件名还是原来的，找不到ige文件，所以就打不开了。但是在使用Erdas和ArcGIS打开该文件时，会正常打开，于是查看GDAL代码，修改部分代码，能够使GDAL正常打开。 修改的代码位置如下，gdal源代码目录/frmts/hfa/hfaband.cpp中 367行处的代码修改为下面的：
使用GDAL创建Erdas格式的金字塔.

在使用Erdas或者ArcGIS打开栅格图像的时候，会创建一个后缀名为rrd的金字塔文件，用于快速显示图像。那么在使用GDAL编写自己的图像算法后，像快速的在Erdas或者ArcGIS中显示，就需要自己创建rrd格式的金字塔文件，这样在打开该图像文件时，打开速度会非常快，在我的电脑上一个2G的img不到一秒钟可以全部加载进来。
查看GDAL中，有个gdaladdo的工具，就是一个专门用于创建金字塔文件的，但是默认的创建的是一个后缀名为ovr的金字塔文件，该种格式的金字塔不能被Erdas或者ArcGIS使用，但是可以在QGIS等使用。为了能够在Erdas中使用，在创建的时候需要指定一个选项，那就是 USE_RRD=YES，使用该选项后，创建的金字塔格式是一个aux文件，虽然后缀名不是rrd，但是该文件是可以被Erdas或者ArcGIS中使用的。
关于gdaladdo的使用帮助，可以参考网址：http://www.gdal.org/gdaladdo.html
Erdas的金字塔是按照2的次方来采样，金字塔顶层的大小应该是小于等于64*64，创建金字塔的，于是按照gdaladdo中的说明，其命令行参数应该是：

gdaladdo --config USE_RRD YES airphoto.img 2 4 8 16 ...

最后的...表示采样级别，一直到最顶层的像元个数小于等于64*64结束。有了上面的知识，下面就给出我写的一个函数，用于创建金字塔。
[cpp] view plaincopyprint?

/**
* @brief 创建金字塔
* @param pszFileName 输入的栅格文件
* @param pProgress 进度条指针
* @return 成功返回0
*/
int CreatePyramids(const char* pszFileName, LT_Progress *pProgress) { if (pProgress != NULL) { pProgress->SetProgressCaption("创建金字塔"); pProgress->SetProgressTip("正在创建金字塔..."); }
GDALAllRegister(); CPLSetConfigOption("USE_RRD","YES"); //创建Erdas格式的字塔文件

/* -------------------------------------------------------------------- */
/* Open data file. */
/* -------------------------------------------------------------------- */

GDALDatasetH hDataset; hDataset = GDALOpen( pszFileName, GA_ReadOnly );
GDALDriverH hDriver = GDALGetDatasetDriver(hDataset); const char* pszDriver = GDALGetDriverShortName(hDriver); if (EQUAL(pszDriver, "HFA") || EQUAL(pszDriver, "PCIDSK")) { GDALClose(hDataset); //如果文件是Erdas的img或者PCI的pix格式，创建内金字塔，其他的创建外金字塔
hDataset = GDALOpen( pszFileName, GA_Update ); }
if( hDataset == NULL ) { if (pProgress != NULL) pProgress->SetProgressTip("打开图像失败，请检查图像是否存在或文件是否是图像文件！");
return RE_NOFILE; }
/* -------------------------------------------------------------------- */
/* Get File basic infomation */
/* -------------------------------------------------------------------- */
int iWidth = GDALGetRasterXSize(hDataset); int iHeigh = GDALGetRasterYSize(hDataset);
int iPixelNum = iWidth * iHeigh; //图像中的总像元个数
int iTopNum = 4096; //顶层金字塔大小，64*64
int iCurNum = iPixelNum / 4;
int anLevels[1024] = { 0 }; int nLevelCount = 0; //金字塔级数

do //计算金字塔级数，从第二级到顶层
{ anLevels[nLevelCount] = static_cast<int>(pow(2.0, nLevelCount+2)); nLevelCount ++; iCurNum /= 4; } while (iCurNum > iTopNum);
const char *pszResampling = "nearest"; //采样方式
GDALProgressFunc pfnProgress = GDALProgress;//进度条

/* -------------------------------------------------------------------- */
/* Generate overviews. */
/* -------------------------------------------------------------------- */
if (nLevelCount > 0 && GDALBuildOverviews( hDataset,pszResampling, nLevelCount, anLevels, 0, NULL, pfnProgress, pProgress ) != CE_None ) { if (pProgress != NULL) pProgress->SetProgressTip("创建金字塔失败！");
return RE_FAILD; }
/* -------------------------------------------------------------------- */
/* Cleanup */
/* -------------------------------------------------------------------- */
GDALClose(hDataset); GDALDestroyDriverManager();
if (pProgress != NULL) pProgress->SetProgressTip("创建金字塔完成！");
return RE_SUCCESS; }

需要说明的是，这段代码创建img格式和pix格式的金字塔会创建内金字塔，Erdas的img格式和PCI的pix格式可以把金字塔存放在文件内部。
PS：在给img文件创建内金字塔后，使用除ArcGIS9.2以外的软件打开后，都正常，但是使用ArcGIS9.2打开后会出现图层偏移的问题，不知道是否ArcGIS9.2的bug。ArcGIS10正常！ArcGIS9.3没有测试。
测试代码：
[cpp] view plaincopyprint?

void main() { LT_ConsoleProgress *pProgress = new LT_ConsoleProgress();
int f = CreatePyramids("C://Work//Data//ttttt.img", pProgress);
if (f == RE_SUCCESS) printf("计算成功/n"); else
printf("计算失败/n");
delete pProgress; }

测试：在使用ArcGIS10打开没有金字塔的文件时，提示：



运行测试代码：



再次打开刚才的文件，没有上面的提示对话框了，而且很快加载进来，说明已经有金字塔了，如果使用Erdas打开的话，可以看到详细的金字塔信息。不过可以试用QGIS查看金字塔信息。右侧的列表显示的是金字塔的级别，Erdas的金字塔是从第二级开始建立的，所以看到第一级的图标上有个红色的小叉叉。见下图：




如何使用GDAL重采样图像
在编写重采样图像时，可以使用GDAL来读写图像，然后自己编写重采样算法（最邻近像元法，双线性内插法，三次立方卷积法等）【关于这采样算法有时间我会单独写一篇文章来说明原理的】将计算的结果写入图像中来实现。
在GDAL的算法中，已经提供了五种重采样算法，其定义如下（位置gdalwarper.h 的46行）：
[cpp] view plaincopyprint?

/*! Warp Resampling Algorithm */
typedef enum { /*! Nearest neighbour (select on one input pixel) */ GRA_NearestNeighbour=0, /*! Bilinear (2x2 kernel) */ GRA_Bilinear=1, /*! Cubic Convolution Approximation (4x4 kernel) */ GRA_Cubic=2, /*! Cubic B-Spline Approximation (4x4 kernel) */ GRA_CubicSpline=3, /*! Lanczos windowed sinc interpolation (6x6 kernel) */ GRA_Lanczos=4 } GDALResampleAlg;

在查看Gdalwarp的源代码发现，warp的功能非常强大，可以用来做投影转换，重投影，投影定义，重采样，镶嵌，几何精校正和影像配准等。一句话，很好很强大。下面就看看其中的一点点皮毛，使用warp来编写一个重采样的接口，代码如下：
[cpp] view plaincopyprint?

/**
* 重采样函数(GDAL)
* @param pszSrcFile 输入文件的路径
* @param pszOutFile 写入的结果图像的路径
* @param fResX X转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小
* @param fResY Y转换采样比，默认大小为1.0
* @param nResampleMode 采样模式，有五种，具体参见GDALResampleAlg定义，默认为双线性内插
* @param pExtent 采样范围，为NULL表示计算全图
* @param pBandIndex 指定的采样波段序号，为NULL表示采样全部波段
* @param pBandCount 采样的波段个数，同pBandIndex一同使用，表示采样波段的个数
* @param pszFormat 写入的结果图像的格式
* @param pProgress 进度条指针
* @return 成功返回0，否则为其他值
*/
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, float fResX , float fResY, LT_ResampleMode nResampleMode, LT_Envelope* pExtent, int* pBandIndex, int *pBandCount, const char* pszFormat, LT_Progress *pProgress) { if(pProgress != NULL) { pProgress->SetProgressCaption("重?采?样?"); pProgress->SetProgressTip("正?在?执?行?重?采?样?..."); }
GDALAllRegister(); GDALDataset *pDSrc = (GDALDataset *)GDALOpen(pszSrcFile, GA_ReadOnly); if (pDSrc == NULL) { if(pProgress != NULL) pProgress->SetProgressTip("指?定?的?文?件?不?存?在?，?或?者?该?格?式?不?被?支?持?！?");
return RE_NOFILE; }
GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); if (pDriver == NULL) { if(pProgress != NULL) pProgress->SetProgressTip("不?能?创?建?该?格?式?的?文?件?！?");
GDALClose((GDALDatasetH) pDSrc); return RE_CREATEFILE; }
int iBandCount = pDSrc->GetRasterCount(); string strWkt = pDSrc->GetProjectionRef(); GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
double dGeoTrans[6] = {0}; pDSrc->GetGeoTransform(dGeoTrans);
int iNewBandCount = iBandCount; if (pBandIndex != NULL && pBandCount != NULL) { int iMaxBandIndex = pBandIndex[0]; //找?出?最?大?的?波?段?索?引?序?号?
for (int i=1; i<*pBandCount; i++) { if (iMaxBandIndex < pBandIndex[i]) iMaxBandIndex = pBandIndex[i]; }
if(iMaxBandIndex > iBandCount) { if(pProgress != NULL) pProgress->SetProgressTip("指?定?的?波?段?序?号?超?过?图?像?的?波?段?数?，?请?检?查?输?入?参?数?！?");
GDALClose((GDALDatasetH) pDSrc); return RE_PARAMERROR; }
iNewBandCount = *pBandCount; }
LT_Envelope enExtent; enExtent.setToNull();
if (pExtent == NULL) //全?图?计?算?
{ double dPrj[4] = {0}; //x1,x2,y1,y2
ImageRowCol2Projection(dGeoTrans, 0, 0, dPrj[0], dPrj[2]); ImageRowCol2Projection(dGeoTrans, pDSrc->GetRasterXSize(), pDSrc->GetRasterYSize(), dPrj[1], dPrj[3]); enExtent.init(dPrj[0], dPrj[1], dPrj[2], dPrj[3]);
pExtent = &enExtent; }
dGeoTrans[0] = pExtent->getMinX(); dGeoTrans[3] = pExtent->getMaxY(); dGeoTrans[1] = dGeoTrans[1] / fResX; dGeoTrans[5] = dGeoTrans[5] / fResY;
int iNewWidth = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[1]) + 0.5) ); int iNewHeight = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[5]) + 0.5) );
GDALDataset *pDDst = pDriver->Create(pszOutFile, iNewWidth, iNewHeight, iNewBandCount, dataType, NULL); if (pDDst == NULL) { if(pProgress != NULL) pProgress->SetProgressTip("创?建?输?出?文?件?失?败?！?");
GDALClose((GDALDatasetH) pDSrc); return RE_CREATEFILE; }
pDDst->SetProjection(strWkt.c_str()); pDDst->SetGeoTransform(dGeoTrans);
GDALResampleAlg eResample = (GDALResampleAlg) nResampleMode;
if(pProgress != NULL) { pProgress->SetProgressTip("正?在?执?行?重?采?样?..."); pProgress->SetProgressTotalStep(iNewBandCount*iNewHeight); }
int *pSrcBand = NULL; int *pDstBand = NULL; int iBandSize = 0; if (pBandIndex != NULL && pBandCount != NULL) { iBandSize = *pBandCount; pSrcBand = new int[iBandSize]; pDstBand = new int[iBandSize];
for (int i=0; i<iBandSize; i++) { pSrcBand[i] = pBandIndex[i]; pDstBand[i] = i+1; } } else
{ iBandSize = iBandCount; pSrcBand = new int[iBandSize]; pDstBand = new int[iBandSize];
for (int i=0; i<iBandSize; i++) { pSrcBand[i] = i+1; pDstBand[i] = i+1; } }
void *hTransformArg = NULL, *hGenImgPrjArg = NULL; hTransformArg = hGenImgPrjArg = GDALCreateGenImgProjTransformer2((GDALDatasetH) pDSrc, (GDALDatasetH) pDDst, NULL); if (hTransformArg == NULL) { if(pProgress != NULL) pProgress->SetProgressTip("转?换?参?数?错?误?！?");
GDALClose((GDALDatasetH) pDSrc); GDALClose((GDALDatasetH) pDDst); return RE_PARAMERROR; }
GDALTransformerFunc pFnTransformer = GDALGenImgProjTransform; GDALWarpOptions *psWo = GDALCreateWarpOptions();
psWo->papszWarpOptions = CSLDuplicate(NULL); psWo->eWorkingDataType = dataType; psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH) pDSrc; psWo->hDstDS = (GDALDatasetH) pDDst;
psWo->pfnTransformer = pFnTransformer; psWo->pTransformerArg = hTransformArg;
psWo->pfnProgress = GDALProgress; psWo->pProgressArg = pProgress;
psWo->nBandCount = iNewBandCount; psWo->panSrcBands = (int *) CPLMalloc(iNewBandCount*sizeof(int)); psWo->panDstBands = (int *) CPLMalloc(iNewBandCount*sizeof(int)); for (int i=0; i<iNewBandCount; i++) { psWo->panSrcBands[i] = pSrcBand[i]; psWo->panDstBands[i] = pDstBand[i]; }
RELEASE(pSrcBand); RELEASE(pDstBand);
GDALWarpOperation oWo; if (oWo.Initialize(psWo) != CE_None) { if(pProgress != NULL) pProgress->SetProgressTip("转?换?参?数?错?误?！?");
GDALClose((GDALDatasetH) pDSrc); GDALClose((GDALDatasetH) pDDst);
return RE_PARAMERROR; }
oWo.ChunkAndWarpImage(0, 0, iNewWidth, iNewHeight);
GDALDestroyGenImgProjTransformer(psWo->pTransformerArg); GDALDestroyWarpOptions( psWo ); GDALClose((GDALDatasetH) pDSrc); GDALClose((GDALDatasetH) pDDst);
if(pProgress != NULL) pProgress->SetProgressTip("重?采?样?完?成?！?");
return RE_SUCCESS; }