关于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行处的代码修改为下面的：

/* -------------------------------------------------------------------- */
/*      Open raw data file.                                             */
/* -------------------------------------------------------------------- */
    const char *pszRawFilename = poDMS->GetStringField( "fileName.string" );
    const char *pszFullFilename;

    pszFullFilename = CPLFormFilename( psInfo->pszPath, pszRawFilename, NULL );

    if( psInfo->eAccess == HFA_ReadOnly )
        fpExternal = VSIFOpenL( pszFullFilename, "rb" );
    else
        fpExternal = VSIFOpenL( pszFullFilename, "r+b" );

    if( fpExternal == NULL )
    {
        CPLString strFileName = psInfo->pszFilename;
        strFileName = strFileName.substr(strFileName.find_last_of('.')+1) + "ige";
        pszFullFilename = CPLFormFilename( psInfo->pszPath, strFileName.c_str(), NULL );

        if( psInfo->eAccess == HFA_ReadOnly )
            fpExternal = VSIFOpenL( pszFullFilename, "rb" );
        else
            fpExternal = VSIFOpenL( pszFullFilename, "r+b" );

        if( fpExternal == NULL )
        {
            CPLError( CE_Failure, CPLE_OpenFailed, 
                      "Unable to open external data file:/n%s/n", 
                      pszFullFilename );
            return CE_Failure;
        }

        psInfo->pszIGEFilename = const_cast<char*>(strFileName.c_str());
    }

希望对大家有用！

使用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 plain copy print ?

/**
* @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 plain copy print ?

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 plain copy print ?

/*! 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 plain copy print ?

/**
* 重采样函数(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;
}