最近有个项目,其中有个功能是要将遥感影像按标准图幅分割,一开始用AE的接口,慢的让人抓狂,就改用GDAL,速度提升很大。我主要通过http://blog.csdn.net/liminlu0314/学习GDAL。本篇主要记录GDAL实现分割的代码,下篇用AE写个demo。
int CutImageByGDAL(const char* pszInFile,const char* pszOutFile,double XMin,double XMax,double YMin,double YMax,const char* pszFormat="GTiff")
{
int iSrcXOffset=,iSrcYOffset=;
int iDstXOffset=,iDstYOffset=;
int ivXSize=,ivYSize=;
int xSize=,ySize=; //结果影像的大小
double fSrcXMin=XMin,fSrcYMax=YMax; GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");
GDALDataset * pSrcDS = (GDALDataset*) GDALOpen(pszInFile, GA_ReadOnly);
if (pSrcDS == NULL)
{
return -;
//return "未能打开影像文件";
}
try
{
int iBandCount = pSrcDS->GetRasterCount();
const char* pszWkt = pSrcDS->GetProjectionRef(); //影像的坐标系
GDALDataType dT = pSrcDS->GetRasterBand()->GetRasterDataType(); //影像的数据类型 double adfGeoTransform[] = {};
double newGeoTransform[] = {};
//获取影像的坐标转换参数
pSrcDS->GetGeoTransform(adfGeoTransform);
//结果影像的坐标转换参数
memcpy(newGeoTransform, adfGeoTransform, sizeof(double)*);
newGeoTransform[]=XMin;
newGeoTransform[]=YMax; //地理坐标转换为影像上的像素坐标
Projection2ImageRowCol(adfGeoTransform,XMin,YMax,iSrcXOffset,iSrcYOffset);
Projection2ImageRowCol(adfGeoTransform,XMax,YMin,xSize,ySize);
xSize=xSize-iSrcXOffset;
ySize=ySize-iSrcYOffset;
ivXSize=xSize;
ivYSize=ySize; GDALDataset* pDstDS;
GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if (pDriver == NULL)
{
GDALClose((GDALDatasetH) pSrcDS);
return -;
//return "未能创建影像文件驱动";
}
//GDALDataset*
pDstDS = pDriver->Create(pszOutFile, xSize, ySize, iBandCount, dT, NULL);
if (pDstDS == NULL)
{
GDALClose((GDALDatasetH) pSrcDS);
return -;
//return "未能创建影像文件";
}
pDstDS->SetGeoTransform(newGeoTransform);
//设置影像坐标系
pDstDS->SetProjection(pszWkt); //边界处理
//判断切割范围,使其不超过原始影像范围
//如果切割范围超过原始影像,会导致结果影像像素值全是NoData
if(XMin<adfGeoTransform[])
fSrcXMin=adfGeoTransform[];
if(YMax>adfGeoTransform[])
fSrcYMax=adfGeoTransform[]; Projection2ImageRowCol(adfGeoTransform,fSrcXMin,fSrcYMax,iSrcXOffset,iSrcYOffset);
Projection2ImageRowCol(newGeoTransform,fSrcXMin,fSrcYMax,iDstXOffset,iDstYOffset); if(iSrcXOffset+xSize>pSrcDS->GetRasterXSize())
ivXSize=pSrcDS->GetRasterXSize()-iSrcXOffset;
if(iDstXOffset+ivXSize>pDstDS->GetRasterXSize())
ivXSize=pDstDS->GetRasterXSize()-iDstXOffset; if(iSrcYOffset+ySize>pSrcDS->GetRasterYSize())
ivYSize=pSrcDS->GetRasterYSize()-iSrcYOffset;
if(iDstYOffset+ivYSize>pDstDS->GetRasterYSize())
ivYSize=pDstDS->GetRasterYSize()-iDstYOffset; void * pMemData;
switch(dT)
{
case GDT_Byte:
pMemData=new char[xSize*ySize*iBandCount];
break;
case GDT_UInt16:
case GDT_Int16:
case GDT_CInt16:
pMemData=new int[xSize*ySize*iBandCount];
break;
case GDT_UInt32:
case GDT_Int32:
case GDT_Float32:
case GDT_CInt32:
case GDT_CFloat32:
pMemData=new float[xSize*ySize*iBandCount];
break;
case GDT_Unknown:
case GDT_Float64:
case GDT_CFloat64:
pMemData=new double[xSize*ySize*iBandCount];
break;
}
pSrcDS->RasterIO(GF_Read,iSrcXOffset,iSrcYOffset,ivXSize,ivYSize,pMemData,ivXSize,ivYSize,dT,iBandCount,,,,);
pDstDS->RasterIO(GF_Write,iDstXOffset,iDstYOffset,ivXSize,ivYSize,pMemData,ivXSize,ivYSize,dT,iBandCount,,,,); GDALClose((GDALDatasetH) pSrcDS);
GDALClose((GDALDatasetH) pDstDS);
free(pMemData); return ;
//return "完成";
}
catch(const char* excep)
{
if(pSrcDS!=NULL)
GDALClose((GDALDatasetH) pSrcDS);
return ;
}
}
Projection2ImageRowCol是把地理坐标转换为影像像素坐标的函数,实现如下
bool Projection2ImageRowCol(double *adfGeoTransform, double dProjX, double dProjY, int &iCol, int &iRow)
{
try
{
double dTemp = adfGeoTransform[]*adfGeoTransform[] - adfGeoTransform[]*adfGeoTransform[];
double dCol = 0.0, dRow = 0.0;
dCol = (adfGeoTransform[]*(dProjX - adfGeoTransform[]) -
adfGeoTransform[]*(dProjY - adfGeoTransform[])) / dTemp +0.5;
dRow = (adfGeoTransform[]*(dProjY - adfGeoTransform[]) -
adfGeoTransform[]*(dProjX - adfGeoTransform[])) / dTemp +0.5; iCol = static_cast<int>(dCol);
iRow = static_cast<int>(dRow);
return true;
}
catch(...)
{
return false;
}
}