使用GDAL对遥感影像进行投影转换

1.在进行投影转换前需要做的工作：

首先编译GDAL和proj4库,GDAL中进行投影转换依赖的开源库是proj4，可以到网上下载已经编译好的包含proj4的GDAL，优快云上就有很多资源。接下来将编译好的GDAL和PROJ4的lib和include路径分别加在VS2010的“附加库目录”，和“附加包含目录”中；在连接器->附加依赖项中添加编译好的proj.lib文件和gdal_i.lib。然后将proj.dll拷到相应的debug文件下。

环境设置好之后，接下来进行投影转换的第一步。

2.投影转换：

我做的是同一椭球体基准面下的转换，将投影坐标转化为地理坐标。只要把原理和环境弄清楚了，相信其他的转换也很简单。

OGRSpatialReference Raster_spf;

OGRSpatialReference *Dest_Raster_spf;

Raster_spf=m_pDataset->GetProjectionRef();  
  if (Raster_spf.IsProjected()) 
  {
    char*pszSrcWKT=NULL;
   Raster_spf.exportToWkt(&pszSrcWKT);
   Dest_Raster_spf=Raster_spf.CloneGeogCS();
    char *pszDestWKT=NULL;
   Dest_Raster_spf->exportToPrettyWkt(&pszDestWKT);
   OGRCoordinateTransformation  *poCT;
   poCT=OGRCreateCoordinateTransformation(&Raster_spf,Dest_Raster_spf);//反过来就是地理坐标转投影坐标
   if (NULL==poCT)
   {
    QMessageBox::warning(this, tr("Failed"), tr("poCT is null!"), tr("OK")); 
   }
   int ndbFlag=poCT->Transform(4,dbX,dbY,NULL);
   if (ndbFlag)
   {
    OGRCoordinateTransformation::DestroyCT(poCT);
   }

其中dbX和dbY存储的是影像四个顶点坐标。

3,转换后处理

转换完之后需要重新计算影像的地理六参数和行列号，然后再利用GDAL中gdalwarp这个接口进行重采样。这部分的内容是借鉴李民录的GDAL书籍和zhouxuguang236的优快云博客，地址为http://blog.youkuaiyun.com/zhouxuguang236/article/details/17468171。

首先，重新计算地理六参数和行列号

dbRes = dbRes / 111000;  //投影坐标转地理坐标分辨率需除以111000m,反过来就得乘以111000m.原理参考上述博客
    double dbMinx = 0; 
    double dbMaxx = 0; 
    double dbMiny = 0; 
    double dbMaxy = 0; 
    dbMinx = min(min(min(dbX[0],dbX[1]),dbX[2]),dbX[3]); 
    dbMaxx = max(max(max(dbX[0],dbX[1]),dbX[2]),dbX[3]); 
    dbMiny = min(min(min(dbY[0],dbY[1]),dbY[2]),dbY[3]); 
    dbMaxy = max(max(max(dbY[0],dbY[1]),dbY[2]),dbY[3]); 
    //估算行列号 
    padfTransform[0] = dbMinx;     //左上角点坐标 
    padfTransform[3] = dbMaxy; 
    //padfTransform[1] 像素宽度, padfTransform[5]像素高度 
    padfTransform[1] = dbRes; 
    padfTransform[5] = -dbRes; 
    //估算行列数 
    int nPixels = ceil(fabs(dbMaxx-dbMinx)/dbRes); 
    int nLines = ceil(fabs(dbMaxy-dbMiny)/dbRes);   

接下来，利用gdalwarp来重采样并将其写入新的影像中。

GDALDriver* pGDalDriver = NULL;//用于creat的驱动
    GDALDataset * poDataset;            //结果影像
    GDALResampleAlg eResampleMethod=GRA_Bilinear;  //默认设置采样方式为双线性插值，时间会稍短一些
    int iDstWidth=nPixels;
    int iDstHeight=nLines; 
    const char *pszFormat;
    pszFormat= m_pDataset->GetDriver()->GetDescription();
    pGDalDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    QString FilenameD;                    
    FilenameD=QFileDialog::getSaveFileName(this,tr("Save Image"),"",tr("Images (*.tif)"));
    if(FilenameD.isEmpty()) 
    {   return; }  
     poDataset=pGDalDriver->Create(FilenameD.toStdString().c_str(),iDstWidth,iDstHeight,dataBands,dataType,NULL);
     poDataset->SetGeoTransform(padfTransform);
     poDataset->SetProjection(pszDestWKT); 
   
     //构造坐标转换关系
     void *hTransformArg=NULL;
     hTransformArg=GDALCreateGenImgProjTransformer((GDALDatasetH)m_pDataset,pszSrcWKT,(GDALDatasetH)poDataset,
      pszDestWKT,FALSE,0.0,1);
     GDALTransformerFunc pfnTransformer=GDALGenImgProjTransform;

     //构造gdalwarp的变化选项
     GDALWarpOptions *psWO=GDALCreateWarpOptions();
     psWO->papszWarpOptions=CSLDuplicate(NULL);
     psWO->eWorkingDataType=dataType;
     psWO->eResampleAlg=eResampleMethod;

     psWO->hSrcDS=(GDALDatasetH)m_pDataset;
     psWO->hDstDS=(GDALDatasetH)poDataset;
     psWO->pfnTransformer=pfnTransformer;
     psWO->pTransformerArg=hTransformArg;

     psWO->nBandCount=dataBands;
     psWO->panSrcBands=(int*)CPLMalloc(psWO->nBandCount*sizeof(int));
     psWO->panDstBands=(int*)CPLMalloc(psWO->nBandCount*sizeof(int));

     psWO->panSrcBands[0]=1;
     psWO->panDstBands[0]=1;

     //创建GDALWarp执行对象并使用GDALWarpOptions来进行初始化
     GDALWarpOperation oWo;
     oWo.Initialize(psWO);

     //执行处理,转换会耗点时间，请耐心等待
     oWo.ChunkAndWarpImage(0,0,iDstWidth,iDstHeight);

     //释放和关闭文件
     GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
     GDALDestroyWarpOptions(psWO);
     GDALClose((GDALDatasetH)m_pDataset);
     GDALClose((GDALDatasetH)poDataset);

  到此，投影转换的整个过程就结束了。之前遇到的问题是在转换过程中执行poCT=OGRCreateCoordinateTransformation(&Raster_spf,Dest_Raster_spf);得到的poCT为NULL，这肯定就是跟proj4库有关了，后来才发现是环境设置的问题。还有在调试时，执行到oWo.ChunkAndWarpImage(0,0,iDstWidth,iDstHeight);就不动了，原以为是转换除了问题，但跟踪发现结果影像中相关的信息都得到了，应该不会有错。后来发现是重采样过程比较慢，需要时间，但最邻近和双线性插值可能需要的时间稍短些。