使用GDAL對HDF資料進行geoloc校正

在上一篇部落格中，大概說了下怎麼使用gdal提供的gdalwarp工具來進行校正處理。其實質與envi的glt校正應該是一樣的。我把gdalwarp的代碼封裝了一下，寫了一個類來進行geoloc處理。希望對大家有用。

先是標頭檔，函數的注釋很詳細，就不多說了。後面的源檔案就是從gdalwarp.cpp中摘錄出來的，有興趣的可以看gdalwarp.cpp ，下面的代碼只是把這個檔案中沒有用到的代碼刪除了，同時有些參數直接寫死了。不太清楚的可以直接看原來的代碼。

/***************************************************************************** Time: 2013-03-02* Project: 遙感處理演算法 * Purpose: 地理尋找表校正* Author:  李民錄* Copyright (c) 2013, liminlu0314@gmail.com* Describe: 地理尋找表校正，用於HDF資料的校正演算法*****************************************************************************/#ifndef IMAGEGEOLOCWARP_H  #define IMAGEGEOLOCWARP_H#include "ImgCore.h"/*** \file ImageGeoLocWarp.h * @brief 地理尋找表校正** 地理尋找表校正，主要用於HDF資料的校正。比如海洋衛星資料等，MODIS資料等。* 一般在該類型的資料中，有兩個Latitude和Longitude的32位浮點數的資料集，* 用來儲存資料的經緯度座標，在校正過程中使用這兩個進行處理*/class IMGALG_API CImageGeoLocWarp{public:/*** @brief 建構函式，地理尋找表校正* @param pszSrcFile未經處理資料路徑* @param pszDstFile結果資料路徑* @param pszFormat結果資料格式* @param pProcess進度條指標*/CImageGeoLocWarp(const char* pszSrcFile, const char* pszDstFile,const char* pszFormat = "HFA", CProcessBase* pProcess = NULL);/*** @brief 解構函式*/~CImageGeoLocWarp();/*** @brief 擷取輸出參數資訊* @param iHeight輸出映像高度* @param iWidth輸出映像寬度* @param padfGeoTransform輸出GeoTransform六參數* @param padfExtent輸出映像四至範圍* @return是否計算成功，以及錯誤碼*/int GetSuggestedWarpOutput(int &iHeight, int &iWidth, double *padfGeoTransform = NULL,double *padfExtent = NULL);/*** @brief 執行校正演算法* @param dResX輸出映像橫向解析度，預設為0，表示自動計算輸出象元大小* @param dResY輸出映像縱向解析度，預設為0，表示自動計算輸出象元大小* @param resampling重採樣方式，預設為雙線性插補* @return是否計算成功，以及錯誤碼*/int DoGeoLocWarp(double dResX = 0.0, double dResY = 0.0, ResampleMethod resampling = RM_Bilinear);private://輸入參數const char* m_pszSrcFile;/*! 未經處理資料路徑 */const char* m_pszDstFile;/*! 結果資料路徑 */const char* m_pszFormat;/*! 結果資料格式 */CProcessBase* m_pProcess;/*! 進度條指標 */ResampleMethod m_Resampling; /*! 重採樣方式 */string m_strWkt;/*! 輸出檔案投影串，這裡只能是WGS84經緯度 */};#endif //IMAGEGEOLOCWARP_H

下面是源檔案，主要是對gdalwarp.cpp進行了精簡，去除了和geoloc無關的代碼。

/***************************************************************************** Time: 2013-03-02* Project: 遙感處理演算法* Purpose: 地理尋找表校正* Author:  李民錄* Copyright (c) 2013, liminlu0314@gmail.com* Describe: 地理尋找表校正，用於HDF資料的校正演算法*****************************************************************************/#include "ImageGeoLocWarp.h"//上面的標頭檔#include "ImageCommon.h"//這個標頭檔就是之前的那個使用gdal對資料進行管理的部落格（用到的函數是刪除臨時檔案）#include "gdal_priv.h"#include "gdalwarper.h"#include "ogr_srs_api.h"CImageGeoLocWarp::CImageGeoLocWarp(const char* pszSrcFile, const char* pszDstFile,const char* pszFormat, CProcessBase* pProcess){m_pszSrcFile = pszSrcFile;/*! 未經處理資料路徑 */m_pszDstFile = pszDstFile;/*! 結果資料路徑 */m_pszFormat = pszFormat;/*! 結果資料格式 */m_pProcess = pProcess;/*! 進度條指標 */m_Resampling = RM_Bilinear; /*! 重採樣方式，預設為雙線性 */m_strWkt = SRS_WKT_WGS84;//投影只能是WGS84}CImageGeoLocWarp::~CImageGeoLocWarp(){}int CImageGeoLocWarp::GetSuggestedWarpOutput(int &iHeight, int &iWidth, double *padfGeoTransform, double *padfExtent){if(m_pszSrcFile == NULL || m_pszDstFile == NULL)return RE_PARAMERROR;GDALAllRegister();GDALDatasetH hSrcDS = GDALOpen( m_pszSrcFile, GA_ReadOnly );//開啟檔案if(hSrcDS == NULL ){if(m_pProcess != NULL)m_pProcess->SetMessage("開啟待糾正影像失敗！");return RE_FILENOTEXIST;}if(m_strWkt.empty())//目標投影不存在{if(m_pProcess != NULL)m_pProcess->SetMessage("目標投影不存在！");GDALClose( hSrcDS );return RE_PARAMERROR;}// 建立座標轉換參數char** papszWarpOptions = NULL;papszWarpOptions = CSLSetNameValue( papszWarpOptions, "METHOD", "GEOLOC_ARRAY" );papszWarpOptions = CSLSetNameValue( papszWarpOptions, "DST_SRS", m_strWkt.c_str() );void *hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszWarpOptions );if(hTransformArg == NULL ){if(m_pProcess != NULL)m_pProcess->SetMessage("座標轉換參數錯誤！");CSLDestroy( papszWarpOptions );GDALClose( hSrcDS );return RE_PARAMERROR;}GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;// 擷取目標檔案大小以及轉換參數資訊double adfDstGeoTransform[6] = {0};double adfExtent[4];int nPixels = 0, nLines = 0;CPLErr eErr = GDALSuggestedWarpOutput2( hSrcDS, psInfo->pfnTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines, adfExtent, 0 );iWidth = nPixels;iHeight = nLines;if(padfGeoTransform != NULL)memcpy(padfGeoTransform, adfDstGeoTransform, sizeof(double)*6);if(padfExtent != NULL)memcpy(padfExtent, adfExtent, sizeof(double)*4);GDALDestroyGenImgProjTransformer(hTransformArg);CSLDestroy( papszWarpOptions );GDALClose( hSrcDS );if(eErr != CE_None)return RE_PARAMERROR;elsereturn RE_SUCCESS;}int CImageGeoLocWarp::DoGeoLocWarp(double dResX, double dResY, ResampleMethod resampling){m_Resampling = resampling;if(m_pszSrcFile == NULL || m_pszDstFile == NULL)return RE_PARAMERROR;if(FileIsUsed(m_pszDstFile, m_pProcess))return RE_FILEISUESED;if(m_pProcess != NULL){m_pProcess->ReSetProcess();m_pProcess->SetMessage("正在執行影像地理校正...");}GDALAllRegister();GDALDatasetH hSrcDS = GDALOpen( m_pszSrcFile, GA_ReadOnly );//開啟檔案if(hSrcDS == NULL ){if(m_pProcess != NULL)m_pProcess->SetMessage("開啟待糾正影像失敗！");return RE_FILENOTEXIST;}if(m_strWkt.empty())//目標投影不存在{if(m_pProcess != NULL)m_pProcess->SetMessage("目標投影不存在！");GDALClose( hSrcDS );return RE_PARAMERROR;}// 建立座標轉換參數char** papszWarpOptions = NULL;papszWarpOptions = CSLSetNameValue( papszWarpOptions, "METHOD", "GEOLOC_ARRAY" );papszWarpOptions = CSLSetNameValue( papszWarpOptions, "DST_SRS", m_strWkt.c_str() );void *hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszWarpOptions );if(hTransformArg == NULL ){if(m_pProcess != NULL)m_pProcess->SetMessage("座標轉換參數錯誤！");CSLDestroy( papszWarpOptions );GDALClose( hSrcDS );return RE_PARAMERROR;}GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;// 擷取目標檔案大小以及轉換參數資訊double adfDstGeoTransform[6] = {0};double adfExtent[4];int nPixels = 0, nLines = 0;CPLErr eErr = GDALSuggestedWarpOutput2( hSrcDS, psInfo->pfnTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines, adfExtent, 0 );if(eErr != CE_None){if(m_pProcess != NULL)m_pProcess->SetMessage("擷取目標檔案大小以及轉換參數資訊錯誤！");GDALDestroyGenImgProjTransformer(hTransformArg);CSLDestroy( papszWarpOptions );GDALClose( hSrcDS );return RE_PARAMERROR;}// 如果使用者指定了輸出映像的解析度if ( dResX != 0.0 || dResY != 0.0 ){// 如果只指定了一個，使用自動計算的結果if ( dResX == 0.0 ) dResX = adfDstGeoTransform[1];if ( dResY == 0.0 ) dResY = adfDstGeoTransform[5];// 確保解析度符號正確if ( dResX < 0.0 ) dResX = -dResX;if ( dResY > 0.0 ) dResY = -dResY;// 計算輸出映像的範圍double minX = adfDstGeoTransform[0];double maxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;double maxY = adfDstGeoTransform[3];double minY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;// 按照使用者指定的解析度來計算映像的輸出大小以及範圍nPixels = ( int )((( maxX - minX ) / dResX ) + 0.5 );nLines  = ( int )((( minY - maxY ) / dResY ) + 0.5 );adfDstGeoTransform[0] = minX;adfDstGeoTransform[3] = maxY;adfDstGeoTransform[1] = dResX;adfDstGeoTransform[5] = dResY;}// 獲得波段數目int nBandCount = GDALGetRasterCount(hSrcDS);// 建立輸出檔案資料類型與源檔案相同GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));// 建立輸出結果檔案GDALDriverH hDriver = GDALGetDriverByName( m_pszFormat);if(hDriver == NULL ){if(m_pProcess != NULL)m_pProcess->SetMessage("不支援該種類型資料！");GDALDestroyGenImgProjTransformer(hTransformArg);CSLDestroy( papszWarpOptions );GDALClose( hSrcDS );return RE_FILENOTSUPPORT;}GDALDatasetH  hDstDS = GDALCreate( hDriver, m_pszDstFile, nPixels, nLines, nBandCount, eDT, NULL );if(hDstDS == NULL){if(m_pProcess != NULL)m_pProcess->SetMessage("建立輸出檔案失敗！");GDALDestroyGenImgProjTransformer(hTransformArg);CSLDestroy( papszWarpOptions );GDALClose( hSrcDS );return RE_CREATEFAILED;}// 設定投影等資訊GDALSetProjection( hDstDS, m_strWkt.c_str() );GDALSetGeoTransform( hDstDS, adfDstGeoTransform );if (hTransformArg != NULL)GDALSetGenImgProjTransformerDstGeoTransform( hTransformArg, adfDstGeoTransform);// 逐個波段擷取顏色表GDALColorTableH hCT;for (int i=1; i<=nBandCount; i++)//band從1開始算{hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS, i) );if( hCT != NULL )GDALSetRasterColorTable( GDALGetRasterBand(hDstDS, i), hCT );}// 建立轉換選項CSLDestroy( papszWarpOptions );//先釋放之前的papszWarpOptions = NULL;papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");GDALTransformerFunc pfnTransformer = NULL;pfnTransformer = GDALGenImgProjTransform;//座標轉換函式GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();psWarpOptions->papszWarpOptions = CSLDuplicate(papszWarpOptions);psWarpOptions->hSrcDS = hSrcDS;psWarpOptions->hDstDS = hDstDS;psWarpOptions->nBandCount = nBandCount;// 申請記憶體psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );  psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );for (int j=0; j<nBandCount; j++){psWarpOptions->panSrcBands[j] = j + 1;  psWarpOptions->panDstBands[j] = j + 1;}psWarpOptions->eResampleAlg = (GDALResampleAlg)resampling;psWarpOptions->pfnProgress = ALGTermProgress;//進度條回呼函數;  psWarpOptions->pProgressArg = m_pProcess;//進度條指標psWarpOptions->papszWarpOptions = papszWarpOptions;psWarpOptions->pfnTransformer = pfnTransformer;psWarpOptions->pTransformerArg = hTransformArg;CPLErr Err = CE_None;GDALWarpOperation oOperation;if( oOperation.Initialize( psWarpOptions ) == CE_None )Err = oOperation.ChunkAndWarpImage( 0, 0, nPixels, nLines);else{CSLDestroy( papszWarpOptions );GDALDestroyGenImgProjTransformer( hTransformArg );GDALClose( hDstDS );GDALClose( hSrcDS );RasterDelete(m_pszDstFile);//刪除結果資料return RE_FAILED;}CSLDestroy( papszWarpOptions );GDALDestroyGenImgProjTransformer( hTransformArg );GDALClose( hDstDS );GDALClose( hSrcDS );if (m_pProcess && !m_pProcess->m_bIsContinue){RasterDelete(m_pszDstFile);//刪除結果資料return RE_USERCANCEL;}if (Err != CE_None)return RE_FAILED;return RE_SUCCESS;}

調用代碼如下：

#include <stdio.h>#include <stdlib.h>#include <boost/progress.hpp>//boost計時函數int main(){int iRev = RE_SUCCESS;CConsoleProcess *pPro = new CConsoleProcess();progress_timer *pTime = new progress_timer();  // 開始計時// 對HDF資料中的一個字資料進行處理，下面分別是子資料集的路徑和輸出檔案路徑const char* pszSrc = "HDF4_SDS:UNKNOWN:\"D:\\Data\\hdf\\H1BCLR101106020418636.L2C.HDF\":12";const char* pszDst = "D:\\Data\\hdf\\H1BCLR101106020418636_12.tif";CImageGeoLocWarp warp(pszSrc, pszDst, "HFA", pPro);iRev = warp.DoGeoLocWarp(0.0, 0.0, RM_NearestNeighbour);if(iRev == RE_SUCCESS)printf("success!\n");elseprintf("failed!\n");delete pPro;delete pTime;::system("pause");return 0;}

最後再說一下，代碼中的傳回值，RE_開頭的去我之前的部落格中去找，就是一些常量，用於標記返回的資訊。還有那個進度條的類，也在我之前的部落格都有說明的。上面的代碼中用到的第三方庫就是BOOST和GDAL，演算法中只有GDAL一個，BOOST庫只在main函數中用來計時用的。其他的在我之前的部落格都有說明的。

不要指望把My Code拿過去直接編譯就能過去，那樣的話，我覺得對你也沒有啥進步。如果一個人看懂My Code的話，那些不認識的函數也就大概能猜出是個啥意思了，自己寫一個也很簡單的。學習別人的代碼，不要指望直接就能用，關鍵是還要看別人的思路，然後結合自己的項目（工程）來進行改造來達到自己的目的，這才是一種對自己有利的學習方式。廢話有點多，覺得說的沒道理的後面這些話就當沒有吧。