重采样包括:上采样和下采样;在使用的过程中直接输入重采样后的长、宽即可!!!
#include<iostream>
#include<ogrsf_frmts.h>
#include<ogr_geometry.h>
#include<gdal_priv.h>
#include<gdal.h>
#include "gdalwarper.h"
#include<opencv2/opencv.hpp>
using namespace cv;
using namespace std;
#define RAM_RIZE 200
//重采样到指定长宽
int reSampleGDAL(const char* pszSrcFile, const char* pszOutFile, int newWidth, int newHeight, GDALResampleAlg eResample = GRA_Bilinear)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_Update);
if (pDSrc == NULL)
{
cout << "打开文件失败!" << endl;
return -1;
}
GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (pDriver == NULL)
{
GDALClose((GDALDatasetH)pDSrc);
cout << "创建文件失败!" << endl;
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; //原影响的宽
int nNewHeight = height; //原影像的高
pDSrc->GetGeoTransform(dGeoTrans);
bool bNoGeoRef = false;
double dOldGeoTrans0 = dGeoTrans[0];
//如果没有投影,人为设置一个
if (strlen(pszSrcWKT) <= 0)
{
OGRSpatialReference oSRS;
oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pszSrcWKT);
pDSrc->SetProjection(pszSrcWKT);
dGeoTrans[0] = 30.0;
dGeoTrans[3] = 30.0;
dGeoTrans[1] = 0.00001;
dGeoTrans[5] = -0.00001;
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 */
/*************** 建立重采样后的尺寸后获取投影信息 **************/
float fResX = (1.0 * newWidth) / width;
float fResY = (1.0 * newHeight) / height;
//获取重采样信息
dGeoTrans[1] = dGeoTrans[1] / fResX;
dGeoTrans[5] = dGeoTrans[5] / fResY;
nNewWidth = static_cast<int>(newWidth + 0.5);
nNewHeight = static_cast<int>(newHeight + 0.5);
/*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);
cout << "创建文件失败!" << endl;
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);
cout << "处理过程中出错!" << endl;
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);
cout << "处理过程中出错!" << endl;
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;
}
void main()
{
double sample = (double)getTickCount();
ResampleGDAL_3("..\\testImage\\scs.tif", ".\\result\\scs2.tif", 4981, 5134, GRA_Cubic);
double sample_end = ((double)getTickCount() - sample) / getTickFrequency();
cout << "降采样使用时间:" << sample_end << endl;
cout << "降采样结束!" << endl;
}