最近有个项目，其中有个功能是要将遥感影像按标准图幅分割，一开始用AE的接口，慢的让人抓狂，就改用GDAL，速度提升很大。我主要通过http://blog.csdn.net/liminlu0314/学习GDAL。本篇主要记录GDAL实现分割的代码，下篇用AE写个demo。
1 int CutImageByGDAL(const char* pszInFile,const char* pszOutFile,double XMin,double XMax,double YMin,double YMax,const char* pszFormat="GTiff") 2 { 3 int iSrcXOffset=0,iSrcYOffset=0; 4 int iDstXOffset=0,iDstYOffset=0; 5 int ivXSize=0,ivYSize=0; 6 int xSize=0,ySize=0; //结果影像的大小 7 double fSrcXMin=XMin,fSrcYMax=YMax; 8 9 GDALAllRegister(); 10 CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO"); 11 GDALDataset * pSrcDS = (GDALDataset*) GDALOpen(pszInFile, GA_ReadOnly); 12 if (pSrcDS == NULL) 13 { 14 return -1; 15 //return "未能打开影像文件"; 16 } 17 try 18 { 19 int iBandCount = pSrcDS->GetRasterCount(); 20 const char* pszWkt = pSrcDS->GetProjectionRef(); //影像的坐标系 21 GDALDataType dT = pSrcDS->GetRasterBand(1)->GetRasterDataType(); //影像的数据类型 22 23 double adfGeoTransform[6] = {0}; 24 double newGeoTransform[6] = {0}; 25 //获取影像的坐标转换参数 26 pSrcDS->GetGeoTransform(adfGeoTransform); 27 //结果影像的坐标转换参数 28 memcpy(newGeoTransform, adfGeoTransform, sizeof(double)*6); 29 newGeoTransform[0]=XMin; 30 newGeoTransform[3]=YMax; 31 32 //地理坐标转换为影像上的像素坐标 33 Projection2ImageRowCol(adfGeoTransform,XMin,YMax,iSrcXOffset,iSrcYOffset); 34 Projection2ImageRowCol(adfGeoTransform,XMax,YMin,xSize,ySize); 35 xSize=xSize-iSrcXOffset; 36 ySize=ySize-iSrcYOffset; 37 ivXSize=xSize; 38 ivYSize=ySize; 39 40 GDALDataset* pDstDS; 41 GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); 42 if (pDriver == NULL) 43 { 44 GDALClose((GDALDatasetH) pSrcDS); 45 return -2; 46 //return "未能创建影像文件驱动"; 47 } 48 //GDALDataset* 49 pDstDS = pDriver->Create(pszOutFile, xSize, ySize, iBandCount, dT, NULL); 50 if (pDstDS == NULL) 51 { 52 GDALClose((GDALDatasetH) pSrcDS); 53 return -3; 54 //return "未能创建影像文件"; 55 } 56 pDstDS->SetGeoTransform(newGeoTransform); 57 //设置影像坐标系 58 pDstDS->SetProjection(pszWkt); 59 60 //边界处理 61 //判断切割范围，使其不超过原始影像范围 62 //如果切割范围超过原始影像，会导致结果影像像素值全是NoData 64 if(XMin<adfGeoTransform[0]) 65 fSrcXMin=adfGeoTransform[0]; 66 if(YMax>adfGeoTransform[3]) 67 fSrcYMax=adfGeoTransform[3]; 68 69 Projection2ImageRowCol(adfGeoTransform,fSrcXMin,fSrcYMax,iSrcXOffset,iSrcYOffset); 70 Projection2ImageRowCol(newGeoTransform,fSrcXMin,fSrcYMax,iDstXOffset,iDstYOffset); 71 72 if(iSrcXOffset+xSize>pSrcDS->GetRasterXSize()) 73 ivXSize=pSrcDS->GetRasterXSize()-iSrcXOffset; 74 if(iDstXOffset+ivXSize>pDstDS->GetRasterXSize()) 75 ivXSize=pDstDS->GetRasterXSize()-iDstXOffset; 76 77 if(iSrcYOffset+ySize>pSrcDS->GetRasterYSize()) 78 ivYSize=pSrcDS->GetRasterYSize()-iSrcYOffset; 79 if(iDstYOffset+ivYSize>pDstDS->GetRasterYSize()) 80 ivYSize=pDstDS->GetRasterYSize()-iDstYOffset; 81 82 void * pMemData; 83 switch(dT) 84 { 85 case GDT_Byte: 86 pMemData=new char[xSize*ySize*iBandCount]; 87 break; 88 case GDT_UInt16: 89 case GDT_Int16: 90 case GDT_CInt16: 91 pMemData=new int[xSize*ySize*iBandCount]; 92 break; 93 case GDT_UInt32: 94 case GDT_Int32: 95 case GDT_Float32: 96 case GDT_CInt32: 97 case GDT_CFloat32: 98 pMemData=new float[xSize*ySize*iBandCount]; 99 break; 100 case GDT_Unknown: 101 case GDT_Float64: 102 case GDT_CFloat64: 103 pMemData=new double[xSize*ySize*iBandCount]; 104 break; 105 } 106 107 pSrcDS->RasterIO(GF_Read,iSrcXOffset,iSrcYOffset,ivXSize,ivYSize,pMemData,ivXSize,ivYSize,dT,iBandCount,0,0,0,0); 108 pDstDS->RasterIO(GF_Write,iDstXOffset,iDstYOffset,ivXSize,ivYSize,pMemData,ivXSize,ivYSize,dT,iBandCount,0,0,0,0); 109 110 GDALClose((GDALDatasetH) pSrcDS); 111 GDALClose((GDALDatasetH) pDstDS); 112 free(pMemData); 113 114 return 1; 115 //return "完成"; 116 } 117 catch(const char* excep) 118 { 119 if(pSrcDS!=NULL) 120 GDALClose((GDALDatasetH) pSrcDS); 121 return 0; 122 } 123 }
Projection2ImageRowCol是把地理坐标转换为影像像素坐标的函数，实现如下
1 bool Projection2ImageRowCol(double *adfGeoTransform, double dProjX, double dProjY, int &iCol, int &iRow) 2 { 3 try 4 { 5 double dTemp = adfGeoTransform[1]*adfGeoTransform[5] - adfGeoTransform[2]*adfGeoTransform[4]; 6 double dCol = 0.0, dRow = 0.0; 7 dCol = (adfGeoTransform[5]*(dProjX - adfGeoTransform[0]) - 8 adfGeoTransform[2]*(dProjY - adfGeoTransform[3])) / dTemp +0.5; 9 dRow = (adfGeoTransform[1]*(dProjY - adfGeoTransform[3]) - 10 adfGeoTransform[4]*(dProjX - adfGeoTransform[0])) / dTemp +0.5; 11 12 iCol = static_cast<int>(dCol); 13 iRow = static_cast<int>(dRow); 14 return true; 15 } 16 catch(...) 17 { 18 return false; 19 } 20 }