简述

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

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

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

C++代码

代码基于当前仓库的master分支构建。

// g++ gdal_pansharpen.cpp -o gdal_pansharpen -I../include -L../lib -lgdal#include 
    
     #include 
     
      #include 
      
       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
       
        (        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
        
         (        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;}

效果对比

GDAL融合效果和原始多光谱波段对比

GDAL融合效果和原始全色波段对比

ARCGIS融合效果与原始全色和多光谱对比

ArcGIS融合过程使用工具箱-->系统工具箱-->Data Management Tools-->栅格-->栅格处理-->创建全色锐化的栅格数据集。

ARCGIS融合效果与原始全色

左边ArcGIS融合效果，右边原始多光谱影像

GDAL融合效果与ArcGIS融合效果对比

左边GDAL融合效果，右边ArcGIS融合效果

左边ArcGIS融合效果，右边GDAL融合效果

转载地址：http://pwwgl.baihongyu.com/

你可能感兴趣的文章