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

  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 }