AE + GDAL实现影像按标准图幅分割（上）

　　最近有个项目，其中有个功能是要将遥感影像按标准图幅分割，一开始用AE的接口，慢的让人抓狂，就改用GDAL，速度提升很大。我主要通过http://blog.csdn.net/liminlu0314/学习GDAL。本篇主要记录GDAL实现分割的代码，下篇用AE写个demo。

　　

intCutImageByGDAL(constchar*pszInFile,constchar*pszOutFile,doubleXMin,doubleXMax,doubleYMin,doubleYMax,constchar*pszFormat="GTiff")
{
intiSrcXOffset=0,iSrcYOffset=0;
intiDstXOffset=0,iDstYOffset=0;
intivXSize=0,ivYSize=0;
intxSize=0,ySize=0;//结果影像的大小
doublefSrcXMin=XMin,fSrcYMax=YMax;

GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");
GDALDataset*pSrcDS=(GDALDataset*)GDALOpen(pszInFile,GA_ReadOnly);
if(pSrcDS==NULL)
{
return-1;
//return"未能打开影像文件";
}
try
{
intiBandCount=pSrcDS->GetRasterCount();
constchar*pszWkt=pSrcDS->GetProjectionRef();//影像的坐标系
GDALDataTypedT=pSrcDS->GetRasterBand(1)->GetRasterDataType();//影像的数据类型

doubleadfGeoTransform[6]={0};
doublenewGeoTransform[6]={0};
//获取影像的坐标转换参数
pSrcDS->GetGeoTransform(adfGeoTransform);
//结果影像的坐标转换参数
memcpy(newGeoTransform,adfGeoTransform,sizeof(double)*6);
newGeoTransform[0]=XMin;
newGeoTransform[3]=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-2;
//return"未能创建影像文件驱动";
}
//GDALDataset*
pDstDS=pDriver->Create(pszOutFile,xSize,ySize,iBandCount,dT,NULL);
if(pDstDS==NULL)
{
GDALClose((GDALDatasetH)pSrcDS);
return-3;
//return"未能创建影像文件";
}
pDstDS->SetGeoTransform(newGeoTransform);
//设置影像坐标系
pDstDS->SetProjection(pszWkt);

//边界处理
//判断切割范围，使其不超过原始影像范围
//如果切割范围超过原始影像，会导致结果影像像素值全是NoData
if(XMin<adfGeoTransform[0])
fSrcXMin=adfGeoTransform[0];
if(YMax>adfGeoTransform[3])
fSrcYMax=adfGeoTransform[3];

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)
{
caseGDT_Byte:
pMemData=newchar[xSize*ySize*iBandCount];
break;
caseGDT_UInt16:
caseGDT_Int16:
caseGDT_CInt16:
pMemData=newint[xSize*ySize*iBandCount];
break;
caseGDT_UInt32:
caseGDT_Int32:
caseGDT_Float32:
caseGDT_CInt32:
caseGDT_CFloat32:
pMemData=newfloat[xSize*ySize*iBandCount];
break;
caseGDT_Unknown:
caseGDT_Float64:
caseGDT_CFloat64:
pMemData=newdouble[xSize*ySize*iBandCount];
break;
}
　　　　
pSrcDS->RasterIO(GF_Read,iSrcXOffset,iSrcYOffset,ivXSize,ivYSize,pMemData,ivXSize,ivYSize,dT,iBandCount,0,0,0,0);
pDstDS->RasterIO(GF_Write,iDstXOffset,iDstYOffset,ivXSize,ivYSize,pMemData,ivXSize,ivYSize,dT,iBandCount,0,0,0,0);

GDALClose((GDALDatasetH)pSrcDS);
GDALClose((GDALDatasetH)pDstDS);
free(pMemData);

return1;
//return"完成";
}
catch(constchar*excep)
{
if(pSrcDS!=NULL)
GDALClose((GDALDatasetH)pSrcDS);
return0;
}
}

　　

Projection2ImageRowCol是把地理坐标转换为影像像素坐标的函数，实现如下

boolProjection2ImageRowCol(double*adfGeoTransform,doubledProjX,doubledProjY,int&iCol,int&iRow)
{
try
{
doubledTemp=adfGeoTransform[1]*adfGeoTransform[5]-adfGeoTransform[2]*adfGeoTransform[4];
doubledCol=0.0,dRow=0.0;
dCol=(adfGeoTransform[5]*(dProjX-adfGeoTransform[0])-
adfGeoTransform[2]*(dProjY-adfGeoTransform[3]))/dTemp+0.5;
dRow=(adfGeoTransform[1]*(dProjY-adfGeoTransform[3])-
adfGeoTransform[4]*(dProjX-adfGeoTransform[0]))/dTemp+0.5;

iCol=static_cast<int>(dCol);
iRow=static_cast<int>(dRow);
returntrue;
}
catch(...)
{
returnfalse;
}
}

ViewCode