GDAL|GDAL影像重采样 GDAL|C++编程

学习过程中，用到重采样，感觉不错，转载过来供大家学习，原文地址：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 fResXX转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小。数值上等于采样后图像的宽度和采样前图像宽度的比 * @param fResYY转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小。数值上等于采样后图像的高度和采样前图像高度的比 * @retrieve0成功 * @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 = https://www.it610.com/article/(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(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(nNewWidth*fResX+0.5); nNewHeight = static_cast(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; ipanSrcBands[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了。
改进： 【GDAL|GDAL影像重采样】

/*** * 遥感影像重采样 * @param pszSrcFile输入文件的路径 * @param pszOutFile写入的结果图像的路径 * @param eResample采样模式，有五种，具体参见GDALResampleAlg定义，默认为双线性内插 GRA_NearestNeighbour=0最近邻法，算法简单并能保持原光谱信息不变；缺点是几何精度差，灰度不连续，边缘会出现锯齿状 GRA_Bilinear=1双线性法，计算简单，图像灰度具有连续性且采样精度比较精确；缺点是会丧失细节； GRA_Cubic=2三次卷积法，计算量大，图像灰度具有连续性且采样精度高； GRA_CubicSpline=3三次样条法，灰度连续性和采样精度最佳； GRA_Lanczos=4分块兰索斯法，由匈牙利数学家、物理学家兰索斯法创立，实验发现效果和双线性接近； * @param fResXX转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小。数值上等于采样后图像的宽度和采样前图像宽度的比 * @param fResYY转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小。数值上等于采样后图像的高度和采样前图像高度的比 * @retrieve0成功 * @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 = https://www.it610.com/article/(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(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(nNewWidth*fResX+0.5); nNewHeight = static_cast(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; ipanSrcBands[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; }

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