1. 程式人生 > >使用GDAL對遙感影像進行投影轉換

使用GDAL對遙感影像進行投影轉換

1.在進行投影轉換前需要做的工作:

首先編譯GDAL和proj4庫,GDAL中進行投影轉換依賴的開源庫是proj4,可以到網上下載已經編譯好的包含proj4的GDAL,CSDN上就有很多資源。接下來將編譯好的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,轉換後處理

首先,重新計算地理六引數和行列號

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);就不動了,原以為是轉換除了問題,但跟蹤發現結果影像中相關的資訊都得到了,應該不會有錯。後來發現是重取樣過程比較慢,需要時間,但最鄰近和雙線性插值可能需要的時間稍短些。