GDAL影像重采样

学习过程中,用到重采样,感觉不错,转载过来供大家学习,原文地址:http://blog.csdn.net/giselite/article/details/7792620

    #include "gdal_priv.h"  
    #include "ogrsf_frmts.h"  
    #include "gdalwarper.h"  
      
    /***  
    * 遥感影像重采样   (要求影像必须有投影,否则走不通) 
    * @param pszSrcFile        输入文件的路径  
    * @param pszOutFile        写入的结果图像的路径  
    * @param eResample         采样模式,有五种,具体参见GDALResampleAlg定义,默认为双线性内插  
    * @param fResX             X转换采样比,默认大小为1.0,大于1图像变大,小于1表示图像缩小。数值上等于采样后图像的宽度和采样前图像宽度的比  
    * @param fResY             Y转换采样比,默认大小为1.0,大于1图像变大,小于1表示图像缩小。数值上等于采样后图像的高度和采样前图像高度的比 
    * @retrieve     0   成功 
    * @retrieve     -1  打开源文件失败 
    * @retrieve     -2  创建新文件失败 
    * @retrieve     -3  处理过程中出错 
    */    
    int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, float fResX = 1.0 ,float fResY  = 1.0,GDALResampleAlg eResample = GRA_Bilinear)    
    {    
        GDALAllRegister();    
        CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");     
        GDALDataset *pDSrc = (GDALDataset *)GDALOpen(pszSrcFile, GA_ReadOnly);    
        if (pDSrc == NULL)    
        {    
            return -1;    
        }    
      
        GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");    
        if (pDriver == NULL)    
        {    
            GDALClose((GDALDatasetH) pDSrc);    
            return -2;    
        }    
        int width=pDSrc->GetRasterXSize();    
        int height=pDSrc->GetRasterYSize();  
        int nBandCount = pDSrc->GetRasterCount();     
        GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();    
      
        char *pszSrcWKT = NULL;  
        pszSrcWKT = const_cast<char *>(pDSrc->GetProjectionRef());  
          
        //如果没有投影,人为设置一个    
        if(strlen(pszSrcWKT)<=0)  
        {  
            OGRSpatialReference oSRS;  
            oSRS.SetUTM(50,true);   //北半球  东经120度  
            oSRS.SetWellKnownGeogCS("WGS84");  
            oSRS.exportToWkt(&pszSrcWKT);  
        }  
      
        void *hTransformArg;  
        hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH) pDSrc,pszSrcWKT,NULL,pszSrcWKT,FALSE,0.0,1);  
      
        //(没有投影的影像到这里就走不通了)  
        if (hTransformArg == NULL)    
        {    
            GDALClose((GDALDatasetH) pDSrc);    
            return -3;    
        }  
      
        double dGeoTrans[6] = {0};    
        int nNewWidth=0,nNewHeight=0;  
        if(GDALSuggestedWarpOutput((GDALDatasetH)pDSrc,GDALGenImgProjTransform,hTransformArg,dGeoTrans,&nNewWidth,&nNewHeight)!=CE_None)  
        {  
            GDALClose((GDALDatasetH) pDSrc);    
            return -3;   
        }  
      
        GDALDestroyGenImgProjTransformer(hTransformArg);  
      
        //adfGeoTransform[0] /* top left x */  
        //adfGeoTransform[1] /* w-e pixel resolution */  
        //adfGeoTransform[2] /* rotation, 0 if image is "north up" */  
        //adfGeoTransform[3] /* top left y */  
        //adfGeoTransform[4] /* rotation, 0 if image is "north up" */  
        //adfGeoTransform[5] /* n-s pixel resolution */  
      
        dGeoTrans[1] = dGeoTrans[1] / fResX;  
        dGeoTrans[5] = dGeoTrans[5] / fResY;  
        nNewWidth = static_cast<int>(nNewWidth*fResX+0.5);  
        nNewHeight = static_cast<int>(nNewHeight*fResY+0.5);  
      
        //创建结果数据集  
        GDALDataset *pDDst = pDriver->Create(pszOutFile, nNewWidth, nNewHeight, nBandCount, dataType, NULL);    
        if (pDDst == NULL)    
        {    
            GDALClose((GDALDatasetH) pDSrc);    
            return -2;    
        }    
          
        pDDst->SetProjection(pszSrcWKT);    
        pDDst->SetGeoTransform(dGeoTrans);       
       
        GDALWarpOptions *psWo = GDALCreateWarpOptions();    
      
        //psWo->papszWarpOptions = CSLDuplicate(NULL);    
        psWo->eWorkingDataType = dataType;    
        psWo->eResampleAlg = eResample;    
      
        psWo->hSrcDS = (GDALDatasetH) pDSrc;    
        psWo->hDstDS = (GDALDatasetH) pDDst;    
      
        psWo->pfnTransformer = GDALGenImgProjTransform;    
        psWo->pTransformerArg = GDALCreateGenImgProjTransformer((GDALDatasetH) pDSrc,pszSrcWKT,(GDALDatasetH) pDDst,pszSrcWKT,FALSE,0.0,1);;    
      
        psWo->nBandCount = nBandCount;    
        psWo->panSrcBands = (int *) CPLMalloc(nBandCount*sizeof(int));    
        psWo->panDstBands = (int *) CPLMalloc(nBandCount*sizeof(int));    
        for (int i=0; i<nBandCount; i++)    
        {    
            psWo->panSrcBands[i] = i+1;    
            psWo->panDstBands[i] = i+1;   
        }     
      
        GDALWarpOperation oWo;    
        if (oWo.Initialize(psWo) != CE_None)    
        {    
            GDALClose((GDALDatasetH) pDSrc);    
            GDALClose((GDALDatasetH) pDDst);    
            return -3;    
        }    
      
        oWo.ChunkAndWarpImage(0, 0, nNewWidth, nNewHeight);    
        GDALFlushCache( pDDst );  
      
        GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);    
        GDALDestroyWarpOptions( psWo );       
        GDALClose((GDALDatasetH) pDSrc);    
        GDALClose((GDALDatasetH) pDDst);    
      
        return 0;    
    }   

用法示例:

    ResampleGDAL("F:\\Work\\\\数据\\spline.tif",  
        "F:\\Work\\数据\\Resample_Lanczos.tif",GRA_Lanczos,10,10);  
      
    ResampleGDAL("F:\\Work\\数据\\test.tif",  
        "F:\\Work\\数据\\Resample3.tif",GRA_Lanczos,3,3);  

缺陷:没有投影就不能重采样,关键点在于:GDALCreateGenImgProjTransformer取到的是null,根据gdal源代码,如果没有投影,dGeoTrans[] = {0,1,0,0,0,1},这个时候GDALCreateGenImgProjTransformer就返回null,因此可以将dGeoTrans[0]或dGeoTrans[3]设置为非0,这样GDALCreateGenImgProjTransformer就不返回null了。
改进:

    /***  
    * 遥感影像重采样 
    * @param pszSrcFile        输入文件的路径  
    * @param pszOutFile        写入的结果图像的路径  
    * @param eResample         采样模式,有五种,具体参见GDALResampleAlg定义,默认为双线性内插  
                                GRA_NearestNeighbour=0      最近邻法,算法简单并能保持原光谱信息不变;缺点是几何精度差,灰度不连续,边缘会出现锯齿状     
                                GRA_Bilinear=1              双线性法,计算简单,图像灰度具有连续性且采样精度比较精确;缺点是会丧失细节; 
                                GRA_Cubic=2                 三次卷积法,计算量大,图像灰度具有连续性且采样精度高; 
                                GRA_CubicSpline=3           三次样条法,灰度连续性和采样精度最佳; 
                                GRA_Lanczos=4               分块兰索斯法,由匈牙利数学家、物理学家兰索斯法创立,实验发现效果和双线性接近; 
    * @param fResX             X转换采样比,默认大小为1.0,大于1图像变大,小于1表示图像缩小。数值上等于采样后图像的宽度和采样前图像宽度的比  
    * @param fResY             Y转换采样比,默认大小为1.0,大于1图像变大,小于1表示图像缩小。数值上等于采样后图像的高度和采样前图像高度的比 
    * @retrieve     0   成功 
    * @retrieve     -1  打开源文件失败 
    * @retrieve     -2  创建新文件失败 
    * @retrieve     -3  处理过程中出错 
    */    
    int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, float fResX = 1.0 ,float fResY  = 1.0,GDALResampleAlg eResample = GRA_Bilinear)    
    {    
        GDALAllRegister();    
        CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");     
        GDALDataset *pDSrc = (GDALDataset *)GDALOpen(pszSrcFile, GA_ReadOnly);    
        if (pDSrc == NULL)    
        {    
            return -1;    
        }    
      
      
        GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");    
        if (pDriver == NULL)    
        {    
            GDALClose((GDALDatasetH) pDSrc);    
            return -2;    
        }    
        int width=pDSrc->GetRasterXSize();    
        int height=pDSrc->GetRasterYSize();  
        int nBandCount = pDSrc->GetRasterCount();     
        GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();    
      
      
        char *pszSrcWKT = NULL;  
        pszSrcWKT = const_cast<char *>(pDSrc->GetProjectionRef());  
      
      
        double dGeoTrans[6] = {0};    
        int nNewWidth=width,nNewHeight=height;  
        pDSrc->GetGeoTransform(dGeoTrans);   
      
      
        bool bNoGeoRef = false;  
        double dOldGeoTrans0 = dGeoTrans[0];  
        //如果没有投影,人为设置一个    
        if(strlen(pszSrcWKT)<=0)  
        {  
            //OGRSpatialReference oSRS;  
            //oSRS.SetUTM(50,true); //北半球  东经120度  
            //oSRS.SetWellKnownGeogCS("WGS84");  
            //oSRS.exportToWkt(&pszSrcWKT);  
            //pDSrc->SetProjection(pszSrcWKT);  
      
      
            //////////////////////////////////////////////////////////////////////////  
            dGeoTrans[0]=1.0;  
            pDSrc->SetGeoTransform(dGeoTrans);  
            //////////////////////////////////////////////////////////////////////////  
      
      
            bNoGeoRef = true;  
        }     
      
      
        //adfGeoTransform[0] /* top left x */  
        //adfGeoTransform[1] /* w-e pixel resolution */  
        //adfGeoTransform[2] /* rotation, 0 if image is "north up" */  
        //adfGeoTransform[3] /* top left y */  
        //adfGeoTransform[4] /* rotation, 0 if image is "north up" */  
        //adfGeoTransform[5] /* n-s pixel resolution */  
      
      
        dGeoTrans[1] = dGeoTrans[1] / fResX;  
        dGeoTrans[5] = dGeoTrans[5] / fResY;  
        nNewWidth = static_cast<int>(nNewWidth*fResX+0.5);  
        nNewHeight = static_cast<int>(nNewHeight*fResY+0.5);  
      
      
        //创建结果数据集  
        GDALDataset *pDDst = pDriver->Create(pszOutFile, nNewWidth, nNewHeight, nBandCount, dataType, NULL);    
        if (pDDst == NULL)    
        {     
            GDALClose((GDALDatasetH) pDSrc);    
            return -2;    
        }    
          
        pDDst->SetProjection(pszSrcWKT);    
        pDDst->SetGeoTransform(dGeoTrans);       
      
      
        void *hTransformArg = NULL;  
        hTransformArg = GDALCreateGenImgProjTransformer2((GDALDatasetH) pDSrc, (GDALDatasetH) pDDst, NULL); //GDALCreateGenImgProjTransformer((GDALDatasetH) pDSrc,pszSrcWKT,(GDALDatasetH) pDDst,pszSrcWKT,FALSE,0.0,1);  
      
      
        if (hTransformArg == NULL)    
        {    
            GDALClose((GDALDatasetH) pDSrc);    
            GDALClose((GDALDatasetH) pDDst);    
            return -3;    
        }  
       
        GDALWarpOptions *psWo = GDALCreateWarpOptions();    
      
      
        psWo->papszWarpOptions = CSLDuplicate(NULL);    
        psWo->eWorkingDataType = dataType;    
        psWo->eResampleAlg = eResample;    
      
      
        psWo->hSrcDS = (GDALDatasetH) pDSrc;    
        psWo->hDstDS = (GDALDatasetH) pDDst;    
      
      
        psWo->pfnTransformer = GDALGenImgProjTransform;    
        psWo->pTransformerArg = hTransformArg;  
      
      
        psWo->nBandCount = nBandCount;    
        psWo->panSrcBands = (int *) CPLMalloc(nBandCount*sizeof(int));    
        psWo->panDstBands = (int *) CPLMalloc(nBandCount*sizeof(int));    
        for (int i=0; i<nBandCount; i++)    
        {    
            psWo->panSrcBands[i] = i+1;    
            psWo->panDstBands[i] = i+1;   
        }     
      
      
        GDALWarpOperation oWo;    
        if (oWo.Initialize(psWo) != CE_None)    
        {    
            GDALClose((GDALDatasetH) pDSrc);    
            GDALClose((GDALDatasetH) pDDst);    
            return -3;    
        }    
      
        oWo.ChunkAndWarpImage(0, 0, nNewWidth, nNewHeight);   
      
      
        GDALDestroyGenImgProjTransformer(hTransformArg);    
        GDALDestroyWarpOptions( psWo );       
        if(bNoGeoRef)   
        {  
            dGeoTrans[0]=dOldGeoTrans0;  
            pDDst->SetGeoTransform(dGeoTrans);     
            //pDDst->SetProjection("");  
        }   
        GDALFlushCache( pDDst );  
        GDALClose((GDALDatasetH) pDSrc);    
        GDALClose((GDALDatasetH) pDDst);      
        return 0;    
    }   

改进后:没有投影也可以过去,但是如果没有投影,输出结果会自动带上一个投影,并且这个投影信息去不掉。

转载自:https://blog.csdn.net/sunj92/article/details/51787856

You may also like...