简述

最近在GDAL的代码中看见了gdalpansharpen.cpp，于是简单的尝试了一下。

融合后的效果比较差，这应该是我对这个算法的使用还不熟悉，还有些地方没有弄清楚。这个代码比较新，是2.1版本才加上的，我在看的时候，代码还有一些小问题，已经在github上提交了issuse了。

融合使用的数据是我在网上找到的高分一号的一景数据，先做了校正，形成全色波段TIFF(单波段)和多光谱波段TIFF(4波段)。

C++代码

代码基于当前https://github.com/OSGeo/gdal仓库的master分支构建。

// g++ gdal_pansharpen.cpp -o gdal_pansharpen -I../include -L../lib -lgdal

#include <gdalpansharpen.h>

#include <cpl_auto_close.h>

#include <cpl_error.h>

int main()

{

    GDALAllRegister();

    // 打开全色波段(高分辨率)文件

    GDALDatasetH hPanDset = GDALOpen("/mnt/data/GF1_PMS2_E116.5_N39.6_20130501_L1A0000127213-PAN2_rpc.tiff",GA_ReadOnly);

    CPL_AUTO_CLOSE_WARP(hPanDset,GDALClose);

    VALIDATE_POINTER1(hPanDset,"Open Pansharpen file failed",1);

    // 打开多光谱(低分辨率)文件

    GDALDatasetH hMssDset = GDALOpen("/mnt/data/GF1_PMS2_E116.5_N39.6_20130501_L1A0000127213-MSS2_rpc.tiff",GA_ReadOnly);

    CPL_AUTO_CLOSE_WARP(hMssDset,GDALClose);

    VALIDATE_POINTER1(hPanDset,"Open Spectral Band file failed",1);

    int nSpectralBands = GDALGetRasterCount(hMssDset);

    GDALPansharpenOptions opts;

    opts.ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;     // 超分辨率贝叶斯法,当前仅支持brovery加权

    opts.eResampleAlg    = GRIORA_Cubic;                 // 重采样至全色光谱波段分辨率的算法

    opts.nBitDepth     = 0;                            // 多光谱波段位深度，0表示默认

    opts.nWeightCount    = nSpectralBands;               // 权重系数数组元素个数(与输入多光谱波段数一致)

    double* pWeightCount= static_cast<double*>(

        CPLMalloc(opts.nWeightCount * sizeof(double))); // 权重系数数组

    CPL_AUTO_CLOSE_WARP(pWeightCount,CPLFree);

    opts.padfWeights = pWeightCount;

    opts.padfWeights[0] = 0.334;    // 蓝

    opts.padfWeights[1] = 0.333;    // 绿

    opts.padfWeights[2] = 0.333;    // 红

    opts.padfWeights[3] = 0.0;        // 近红外

    // 设置全色波段(高分辨率)

    opts.hPanchroBand = GDALGetRasterBand(hPanDset,1);

    // 设置多光谱波段

    opts.nInputSpectralBands = nSpectralBands;

    GDALRasterBandH* pInputSpectralBands = static_cast<GDALRasterBandH*>(

        CPLMalloc(sizeof(GDALRasterBandH) * nSpectralBands));

    CPL_AUTO_CLOSE_WARP(pInputSpectralBands,CPLFree);

    opts.pahInputSpectralBands = pInputSpectralBands;

    // std::generatr

    for(int i=0;i< nSpectralBands;++i) {

        opts.pahInputSpectralBands[i] = GDALGetRasterBand(hMssDset,i+1);

    }

    // 设置需要输出到全色波段分辨率的波段

    opts.nOutPansharpenedBands = 4;

    // 这个数组里面存的是pahInputSpectralBands里波段的索引值

    int panOutPansharpenedBands[4] = { 2, 1, 0, 3}; // 红、绿、蓝、近红外

    opts.panOutPansharpenedBands = panOutPansharpenedBands;

    opts.bHasNoData = FALSE;   // 全色和多光谱波段是否具有无效值(NoData值)

    opts.dfNoData   = 0.0;     // 全色和多光谱波段的无效值,也将作为输出的NoData值

    opts.nThreads   = -1;      // 使用的线程数，-1表示使用CPU线程数

    // 设置多光谱波段与全色波段在像素上的移位(保证地理空间位置对齐)

    // 都是相对于全色波段的0,0像素的像素(全色波段像素大小)偏移

    // 也就是两者的0,0像素的地理空间上的偏移量在全色波段分辨率相应的像素数

    double adfGTPan[6];

    GDALGetGeoTransform(hPanDset,adfGTPan);

    double adfGTMss[6];

    GDALGetGeoTransform(hPanDset,adfGTMss);

    opts.dfMSShiftX = (adfGTPan[0] - adfGTMss[0]) / adfGTPan[1];

    opts.dfMSShiftY = (adfGTPan[3] - adfGTMss[3]) / adfGTPan[5];

    GDALPansharpenOperation operation;

    CPLErr err = operation.Initialize(&opts);

    if(err != CE_None) {

        return -2;

    }

    // 创建输出文件(和全色波段一样大)

    int nOutXSize, nOutYSize;

    nOutXSize = GDALGetRasterBandXSize(opts.hPanchroBand);

    nOutYSize = GDALGetRasterBandYSize(opts.hPanchroBand);

    GDALDataType eBufDataType = GDALGetRasterDataType(opts.pahInputSpectralBands[0]);

    eBufDataType = GDT_Float64;

    GDALDriverH hDriver = GDALGetDriverByName("GTiff");

    CPLStringList CreateOption;

    CreateOption.AddNameValue("TILED", "YES");

    CreateOption.AddNameValue("BIGTIFF", "YES");

    CreateOption.AddNameValue("INTERLEAVE", "BAND");

    CreateOption.AddNameValue("COMPRESS", "LZW");   // 中度压缩

    CreateOption.AddNameValue("ZLEVEL", "6");

    GDALDatasetH hOutDset = GDALCreate(hDriver,

                "/mnt/data/GF1_PMS2_E116.5_N39.6_20130501_L1A0000127213.tif",

                nOutXSize, nOutYSize, nSpectralBands, GDT_UInt16,

                CreateOption);

    CPL_AUTO_CLOSE_WARP(hOutDset,GDALClose);

    VALIDATE_POINTER1(hOutDset,"Create Output file error", -3);

    GDALSetGeoTransform(hOutDset, adfGTPan);

    GDALSetProjection(hOutDset,GDALGetProjectionRef(hPanDset));

     void* pData = CPLMalloc(256*256*GDALGetDataTypeSizeBytes(eBufDataType)*nSpectralBands);

     CPL_AUTO_CLOSE_WARP(pData,CPLFree);

    for(int nYOff = 0; nYOff < nOutYSize; nYOff += 256) {

        for(int nXOff = 0; nXOff < nOutXSize; nXOff += 256) {

            int nXSize = std::min(nOutXSize - nXOff,256);

            int nYSize = std::min(nOutYSize - nYOff,256);

            printf("Process[%6d,%6d,%6d,%6d]\r",nXOff,nYOff,nXOff+nXSize,nYOff+nYSize);

            err = operation.ProcessRegion(nXOff,nYOff,nXSize,nYSize,pData,eBufDataType);

            if(err == CE_Failure) {

                CPLError(err,CPLE_AppDefined,"operation.ProcessRegion");

                return -4;

            }

            int panBandMap[4] = { 1, 2, 3, 4};

            err = GDALDatasetRasterIO(hOutDset, GF_Write,

                        nXOff,nYOff,nXSize,nYSize,

                        pData,nXSize,nYSize,

                        eBufDataType,

                        4,panBandMap,

                        0,0,nXSize*nYSize*GDALGetDataTypeSizeBytes(eBufDataType));

        }

    }

    puts("\nPansharpen finished");

    return 0;

}