GDAL源碼剖析（十二）之GDAL Warp API使用說明

一、簡介

本文原文地址：http://www.gdal.org/warptut.html

GDAL Warp API（在檔案gdalwarper.h中定義）是一個高效的進行映像變換的介面。主要由幾何變換函數（GDALTransformerFunc），多種映像重採樣方式，掩碼操作選項等組成。這個介面可以對很大的映像進行處理。

下面說明樣本讓你如何在程式中使用變換API。首先假定你已經熟悉了GDAL的抽象資料模型，以及GDAL的API。

在程式中，首先要初始化一個GDALWarpOptions 結構體的對象，然後使用GDALWarpOptions 的對象來初始化GDALWarpOperation 的對象，最後通過調用GDALWarpKernel 類裡面的GDALWarpOperation::ChunkAndWarpImage() 函數來完成映像的變換。

二、簡單的影像重投影樣本

首先我們以一個映像重投影的例子來入手，需要注意的是，這裡假設輸出結果檔案已經存在，同時這裡沒有對錯誤資訊進行檢查，僅僅示範的最正常的處理流程。

#include "gdalwarper.h"int main(){    GDALDatasetH  hSrcDS, hDstDS;     // 開啟輸入輸出映像    GDALAllRegister();    hSrcDS = GDALOpen("in.tif", GA_ReadOnly );   hDstDS = GDALOpen("out.tif", GA_Update );     // 建立變換選項    GDALWarpOptions*psWarpOptions = GDALCreateWarpOptions();    psWarpOptions->hSrcDS =hSrcDS;   psWarpOptions->hDstDS =hDstDS;    psWarpOptions->nBandCount = 1;   psWarpOptions->panSrcBands =       (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );   psWarpOptions->panSrcBands[0] = 1;   psWarpOptions->panDstBands =       (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );   psWarpOptions->panDstBands[0] = 1;     psWarpOptions->pfnProgress = GDALTermProgress;       // 建立重投影變換函數   psWarpOptions->pTransformerArg =        GDALCreateGenImgProjTransformer( hSrcDS,                                         GDALGetProjectionRef(hSrcDS),                                        hDstDS,                                         GDALGetProjectionRef(hDstDS),                                         FALSE,0.0, 1 );   psWarpOptions->pfnTransformer = GDALGenImgProjTransform;     // 初始化並且執行變換操作    GDALWarpOperationoOperation;    oOperation.Initialize(psWarpOptions );   oOperation.ChunkAndWarpImage( 0, 0,                                  GDALGetRasterXSize( hDstDS),                                  GDALGetRasterYSize( hDstDS) );     GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg );   GDALDestroyWarpOptions( psWarpOptions );     GDALClose( hDstDS );    GDALClose( hSrcDS );    return 0;}

這個例子中首先開啟已經存在的輸入檔案和輸出檔案（in.tif和out.tif）。使用GDALCreateWarpOptions()函數建立一個GDALWarpOptions結構體對象（結構體中的參數會指定一個比較合理的預設值），然後對這個結構體對象設定輸入輸出檔案的控制代碼（就是檔案指標）和要處理的波段。需要注意的是panSrcBands和panDstBands數組需要在外面動態申請，然後在調用GDALDestroyWarpOptions()函數的時候會自動釋放，所以在後面就不需要我們對這兩個數組進行釋放了。GDALTermProgress是一個簡單的控制台進度資訊函數，用來顯示處理進度。

GDALCreateGenImgProjTransformer()函數是用來建立一個從原始映像到結果映像的投影變換參數。我們假設這兩個映像都有合適的四至範圍和座標系統，不使用GCP點。

一旦GDALWarpOptions結構體建立好了，可以使用這個GDALWarpOptions對象來初始化一個GDALWarpOperation的對象，然後再調用函數GDALWarpOperation::ChunkAndWarpImage()進行轉換。在轉換完成之後，轉換選項，資料集等都需要進行釋放。

通常應該在開啟映像之後進行一系列的檢查，然後在建立投影變換參數是要對返回的參數進行檢查（傳回值為NULL表示失敗），最後還要對建立的變換操作進行檢查。上面的例子為了方便，沒有對這些資訊進行檢查，在我們自己的程式中需要對這些進行檢查。

三、其他變換選項

GDALWarpOptions結構體中包含了很多參數來對變換進行設定。下面對一些比較重要的進行列舉說明：

1.  GDALWarpOptions::dfWarpMemoryLimit：設定GDALWarpOperation在處理映像中使用的最大記憶體數。單位為位元，預設值比較保守（比較小），可以根據自己的記憶體大小來進行調整這個值。增加這個值可以協助提高程式的運行效率，但是需要注意記憶體的大小。這個大小、GDAL的緩衝大小，還有你的應用程式以及系統所需要的記憶體加起來要小於系統的記憶體，否則會導致錯誤。一般來說，比如一個記憶體為256MB的系統，這個值最少設定為64MB比較合理。注意，這個值不包括GDAL讀取資料使用的記憶體。

2.  GDALWarpOpations::eResampleAlg：重採樣方式，可用的值有 GRA_NearestNeighbour （最鄰近採樣，預設值，處理速度最快）、GRA_Bilinear（2x2雙線性插補採樣）和 GRA_Cubic（三次立方卷積採樣）。GRA_NearestNeighbour採樣方式一般用在專題圖或者彩色映像中，其他的重採樣方式效果比較好，尤其是在計算中改變映像解析度的演算法中。

3.  GDALWarpOptions::padfSrcNoDataReal：這個數組（每個波段一個值），可以用來指定輸入映像波段的NODATA值，比像的背景值，對於這樣的值，演算法不會參與運算，直接將這個值複製到結果映像中。

4.  GDALWarpOptions::papszWarpOptions：這個是一個字串列表，用來設定映像變換過程中的一些選項，樣子為NAME=VALUE。更多詳細的內容可以參考 GDALWarpOptions::papszWarpOptions的文檔，裡面含有全部的選項，支援的值包括：

• INIT_DEST=[value]或者INIT_DEST=NO_DATA：這個選項用來強制設定結果映像的初始值（所有的波段），初始值為指定的value，或者NODATA值。NODATA值從參數padfDstNoDataReal或者參數padfDstNoDataImag中擷取。如果這個值沒有設定，那麼將會使用原始映像的NODATA值來覆蓋。
• WRITE_FLUSH=YES/NO：這個選項用來強制設定在處理完每一塊後將資料寫入磁碟中。在某些時候，這個選項可以更加安全的寫入結果資料，但是同時會增加更多的磁碟操作。目前這個預設值為NO。

四、建立輸出檔案

在前面的例子中，結果映像是已經存在的。選擇我們將通過預測輸出檔案的範圍和座標系統來建立新的映像。這個操作不是映像變換的特殊API，這個僅僅是變換的一個API。

#include "gdalwarper.h"#include"ogr_spatialref.h" ...     GDALDriverH hDriver;    GDALDataType eDT;    GDALDatasetH hDstDS;    GDALDatasetH hSrcDS;     // 開啟源檔案    hSrcDS = GDALOpen("in.tif", GA_ReadOnly );    CPLAssert( hSrcDS != NULL );       // 建立輸出映像的資料類型和輸入映像第一個波段類型一致    eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));     // 擷取輸出映像的驅動（GeoTIFF格式）    hDriver = GDALGetDriverByName("GTiff" );    CPLAssert( hDriver != NULL );     // 擷取源映像的座標系統    const char *pszSrcWKT, *pszDstWKT = NULL;    pszSrcWKT = GDALGetProjectionRef( hSrcDS);    CPLAssert( pszSrcWKT != NULL &&strlen(pszSrcWKT) > 0 );     // 設定輸出映像的座標系統為UTM 11 WGS84    OGRSpatialReference oSRS;    oSRS.SetUTM( 11, TRUE );    oSRS.SetWellKnownGeogCS( "WGS84");    oSRS.exportToWkt( &pszDstWKT );     // 建立一個變換參數，從源映像的行列號座標到結果映像的地理座標（沒有         //結果行列）的控制代碼，初始值設定為NULL。    void *hTransformArg;    hTransformArg =        GDALCreateGenImgProjTransformer( hSrcDS,pszSrcWKT, NULL, pszDstWKT,                                         FALSE,0, 1 );    CPLAssert( hTransformArg != NULL );     // 擷取輸出檔案的近似地理範圍和解析度    double adfDstGeoTransform[6];    int nPixels=0, nLines=0;    CPLErr eErr;    eErr = GDALSuggestedWarpOutput( hSrcDS,                                 GDALGenImgProjTransform, hTransformArg,                                 adfDstGeoTransform, &nPixels, &nLines );    CPLAssert( eErr == CE_None );     GDALDestroyGenImgProjTransformer(hTransformArg );     // 建立輸出檔案    hDstDS = GDALCreate(hDriver, "out.tif", nPixels, nLines,                       GDALGetRasterCount(hSrcDS),eDT, NULL );       CPLAssert( hDstDS != NULL );     // 寫入投影    GDALSetProjection( hDstDS,pszDstWKT );    GDALSetGeoTransform( hDstDS,adfDstGeoTransform );     // 複製顏色表，如果有的話    GDALColorTableH hCT;    hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1));    if( hCT != NULL )        GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1),hCT );     ... 變換處理之前做的工作...

這裡需要注意的一些邏輯關係：

• 我們需要建立的輸出座標是使用地理座標表示的，而不是行列號座標。同時輸出的座標是在結果座標系統之下的座標，不是原始座標系統下的。
•  GDALSuggestedWarpOutput()函數傳回值adfDstGeoTransform、nPixels
  和nLines這三個值是描述輸出映像的大小和四至範圍，這個四至範圍包含了源映像的所有的像素點。裡面的解析度是根據原始映像估算出來的，輸出映像的解析度一致是橫向和縱向保持一致，不受輸入映像限制。
•  變化要求輸出檔案格式是可以“隨機”寫入的。一般使用Create()（不能用CreateCopy()函數）函數建立非壓縮的映像格式。如果要建立壓縮的映像格式，或者使用CreateCopy()函數建立的話，系統內部會建立一個臨時檔案，然後使用CreateCopy()函數寫入到結果映像中。
•  變換API僅僅處理映像的象元值。所有的顏色表，地理參考以及其他中繼資料資訊必須使用程式自行設定到結果映像中去，變換是不會設定這些資訊的。

五、效能最佳化

下面幾點可以在使用變換API的時候提高處理效率。

1. 增加變換API使用的記憶體，可以使程式執行的更快。這個參數叫GDALWarpOptions::dfWarpMemoryLimit。理論上，越大的塊對於資料I/O效率越高，並且變換的效率也會越高。然而，變換所需的記憶體和GDAL緩衝應該小於機器的記憶體，最多是記憶體的2/3左右。
2. 增加GDAL的緩衝。這個尤其對於在處理非常大的輸入映像很有用。如果所有的輸入輸出映像的資料頻繁的讀取會極大的降低運行效率。使用函數GDALSetCacheMax()來設定GDAL使用的最大緩衝。
3.  使用近似的變換代替精確的變換（精確的變換是每個象元都會計算一次）。下面的代碼用來說明近似變換的使用方式，近似變換要求指定一個變換的誤差dfErrorThreshold，這個誤差單位是輸出映像的象元個數。

4. hTransformArg =    GDALCreateApproxTransformer(GDALGenImgProjTransform,hGenImgProjArg, dfErrorThreshold );pfnTransformer = GDALApproxTransform;

5. 當寫入一個空的輸出檔案，在GDALWarpOptions::papszWarpOptions 對象中使用INIT_DEST選項來進行設定。這樣可以很大的減少磁碟的IO操作。
6. 輸入和輸出檔案使用分Block Storage的檔案格式。區塊檔案格式可以快速的訪問一塊資料，相比來說，普通的大資料量的順序隱藏檔格式在IO操作中要花費更多的時間。
7. 在一次調用中處理所有的波段。一次處理所有的波段比每次處理一個波段要保險。
8. 使用GDALWarpOperation::ChunkAndWarpMulti()方法來代替GDALWarpOperation::ChunkAndWarpImage()方法。這個使用多個獨立的線程來進行IO操作和變換影像處理，能夠提高CPU和IO的效率。執行這個操作，GDAL需要多線程的支援（從GDAL1.8.0開始，預設在Win32平台、Unix平台是支援的，對於之前的版本，需要在配置中進行修改才行）。
9. 重採樣方式，最鄰近象元執行速度最快，雙線性插補次之，三次立方卷次最慢。不要使用根據複雜的重採樣函數。
10. 避免使用複雜的掩碼選項。比如，最鄰近採樣在處理沒有掩碼的8bit資料要比一般的效率高很多。

六、關於掩碼選項

GDALWarpOptions包含了一些處理複雜的掩碼的能力，比如掩碼的有效性，對輸入和輸出資料的掩碼。這些功能還沒有做足夠的測試。其他每個波段的有效掩碼在使用的時候要小心。