如何使用GDAL重采样图像

在编写重采样图像时，可以使用GDAL来读写图像，然后自己编写重采样算法（最邻近像元法，双线性内插法，三次立方卷积法等）【关于这采样算法有时间我会单独写一篇文章来说明原理的】将计算的结果写入图像中来实现。

在GDAL的算法中，已经提供了五种重采样算法，其定义如下（位置gdalwarper.h 的46行）：

/*! 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来编写一个重采样的接口，代码如下：

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

PS：在使用Windows Live Writer来写博客，使用VSPaste插件粘贴代码的时候，发现会在汉字后面加一个“？”号，我懒得修改了，大家就凑合看看吧！

如何使用GDAL重采样图像相关推荐

python gdal 重采样_Python遥感影像重采样
图像重采样就是从高分辨率遥感影像中提取出低分辨率影像,或者从低分辨率影像中提取高分辨率影像的过程.使用python和gdal实现图像重采样的方法有多种,下面介绍两种. (1) 使用ReadAsA ...
gdal读写图像分块处理
转自赵文原文 gdal读写图像分块处理(精华版) Review: 用gdal,感觉还不如直接用C++底层函数对遥感数据进行处理.因为gdal进行太多封装,如果你仅仅只是Geotif等格式进行处理,IO ...
学习笔记——【python】GetGeoTransform()使用，gdal截取图像，使用GDAL进行影像投影坐标、地理坐标、图上坐标的转换
1. GetGeoTransform()使用.gdal截取图像 GetGeoTransform() GeoTransform[0],左上角横坐标(应该是投影坐标) GeoTransform[2],行旋 ...
基于GDAL库图像读写——涉及(tif/tiff/bmp/jpg/png/gif等)多种格式图像的I/O
说来也是神奇,之所以会发这篇文章,还是因为在CSDN发文章中无法上传tif格式的图片造成的,如下. 转入文件夹下,并不显示tif图片,选择所有文件确认一副tiff图片,仍然提示上传错误.不知道CSDN ...
python dataset[trans_python gdal根据图像坐标生成矢量框（含图像坐标转地理坐标）...
要生成矢量框需要将图像坐标转换为地理坐标或者投影坐标,以下代码是生成了满足条件的1000*1000区域对应的矢量框,关键在于红色字体部分. # -*- coding: utf-8 -*- import ...
GDAL创建图像提示Driver xxx does not support XXX creation option的原因
经常在群里有人问,创建图像的时候为什么老是提示下面的信息. CPLError: Driver GTiff does not support DCAP_CREATE creation option Wa ...
关于GDAL计算图像坐标的几个问题
使用GDAL处理地理图像时,不可避免的会遇到一个问题,图像的地理坐标问题,因为有了这个地理坐标,地理图像才和普通图像有了最本质的区别,那么在使用GDAL时,如果处理与地理坐标相关的信息呢?下面进行简单 ...
GDAL对图像文件格式的转换
我们常常在图像处理过程中遇到不同软件或程序要求输入的图像格式不同(有些程序或软件支持的数据格式不是常用的Tiff,Img等数据格式),因此需要对不同的数据格式相互进行转换. 我这里以GTiff(.ti ...
（一）GDAL计算图像坐标与像素坐标之间的关系
1. 原理使用GDAL处理地理图像时,不可避免的会遇到一个问题,图像的地理坐标问题,因为有了这个地理坐标,地理图像才和普通图像有了最本质的区别,那么在使用GDAL时,如果处理与地理坐标相关的信息呢? ...

如何使用GDAL重采样图像

如何使用GDAL重采样图像相关推荐

最新文章

热门文章