http://blog.csdn.net/mengdong_zy/article/details/19043689

GDAL柵格影象操作

    GDAL是一個操作各種柵格和向量（由ogr這個庫實現）地理資料格式的開源庫。包括讀取、寫入、轉換、處理各種柵格和向量資料格式（有些特定的格式對一些操作如寫入等不支援）。

即使不是進行地理遙感方面的應用研究，GDAL也是一個非常有用的庫，因為它可以支援大量我們常見的影象資料，比如jpg，gif之類的。完整的格式清單可以到此連結檢視http://www.gdal.org/formats_list.html。而且已經有包括GoogleEarth在內的很多軟體都是在使用GDAL作為後臺庫。

本文就以VC為開發平臺介紹GDAL對柵格資料的操作方法。      

Include目錄是開發中需要的標頭檔案，lib中是所需要的lib檔案，在VC8中應當將其存放目錄新增到目錄列表中，選擇選單的“工具-選項-專案和解決方案-VC++目錄”，分別在“包含檔案”和“庫檔案”中將此兩個目錄新增進去。

在專案的屬性頁中，選擇“配置屬性-連結器-輸入”，在“附加依賴項”中新增gdal_i-vc8.lib和gdal_id-vc8.lib兩個使用GDAL中需要的靜態庫檔案，或者在程式中新增以下兩行程式碼也可以。

#pragma comment(lib, "gdal_i-vc8.lib")

#pragma comment(lib, "gdal_id-vc8.lib")

       Bin目錄下的動態連結庫檔案應當放置於程式能夠訪問的位置，比如windows\system32中。

       此外，在程式中需要引入的標頭檔案是gdal_priv.h。

 

現在開始用C++來對影象檔案進行操作。

在開啟檔案之前需要首先註冊所需要的驅動程式，一般來說我們可以預設註冊所有支援的格式驅動，所使用的函式是GDALAllRegister()。然後就是開啟檔案操作。

       這裡要說一個資料集的概念，也就是所謂的Dataset。在GDAL中可以說資料的核心就是Dataset，簡單來說可以將Dataset就理解為影象檔案，比如說一個jpeg格式的檔案就是一個資料集，當然其他一些檔案格式可能在一個資料集中包含多於一個檔案，比如可能除了影象資料檔案外還可能會有一些附加資訊檔案等。

在資料集下最重要組成部分就是所謂的波段band，波段可多可少，比如一個RGB真彩色的影象就有3個波段，分別代表紅色綠色和藍色波段，如果是灰度圖，那可能就只有一個波段，而很多遙感影象可能就會多於3個波段。

除了波段外，資料集中還含有影象相關的座標系投影資訊，後設資料資訊等資料。

檔案的開啟使用的是GDALOpen ( const char *  pszFilename, GDALAccess  eAccess  )，pszFilename是檔案路徑，eAccess 是訪問許可權，可以是GA_ReadOnly只讀，也可以是GA_Update來對檔案進行修改。比如我們以只讀模式開啟一個tif檔案：

GDALDataset *poDataset; //資料集物件指標

GDALAllRegister();//註冊驅動

     poDataset = (GDALDataset *) GDALOpen( "c:\\terra335h_EV_250_Aggr500_RefSB_b0.tif", GA_ReadOnly );

     if( poDataset != NULL /*檢查是否正常開啟檔案*/)

     {

         //do something

}

       delete poDataset; //釋放資源

 

在確認poDataset不是NULL的情況下就可以對影象資料集進行操作了。

       首先來看一下有關這個影象的基本資訊，比如長寬和波段數。

       cout<<"RasterXSize:"<<poDataset->GetRasterXSize()<<endl;//x方向長度

     cout<<"RasterYSize:"<<poDataset->GetRasterYSize()<<endl;//y方向長度

     cout<<"RasterCount:"<<poDataset->GetRasterCount()<<endl;//波段數量

       我的這個檔案輸出結果如下：

RasterXSize:5684

RasterYSize:3655

RasterCount:1

       這是個5684×3655的灰度圖。

      

有關的投影資訊可以由GetProjectionRef函式獲得。

cout<<poDataset->GetProjectionRef()<<endl;

       輸出一個字串，是WKT格式的座標投影資訊，形如：

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.2572235630016,AUT

HORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["deg

ree",0.0174532925199433],AUTHORITY["EPSG","4326"]]

      

有關仿射變換和GCP的資訊可以通過GetGeoTransform，GetGCPCount，GetGCPProjection，GetGCPs獲得。比如取得仿射變換的資訊如下：

double adfGeoTransform[6];//仿射引數

poDataset->GetGeoTransform(adfGeoTransform);//獲取引數

然後就可以使用以下公式來計算仿射座標了：

Xgeo = adfGeoTransform[0] + Xpixel* adfGeoTransform[1] + Yline* adfGeoTransform[2]

Ygeo = adfGeoTransform[3] + Xpixel* adfGeoTransform[4] + Yline* adfGeoTransform[5]

       這裡的座標系是以左上角為(0,0)點，向右向下為正向。

       Xgeo表示轉換後的地理x座標，Ygeo是其相應y座標，Xpixel是影象中畫素的列座標，Yline則是行座標。而且Xpixel和Yline的(0,0)座標是影象左上角畫素的左上角，因此左上角畫素的中心座標是Xpixel=0.5，Yline=0.5。

 

       接下來，我們就開始進入到波段處理。波段的獲取使用GetRasterBand函式，引數是波段序號，從1開始。比如還是上面那個圖，只有一個波段，那要得到這個波段的操作如下：

GDALRasterBand *poBand;

poBand = poDataset->GetRasterBand( 1 );

關於波段的資訊也很多，我們撿幾個主要的看看。

首先也是尺寸，在一個Dataset裡所有波段的尺寸都一樣，也就是上面用Dataset獲取的數值，就不多數了。

 

然後是資料型別。不同的影象格式的儲存的資料型別是不一樣的，比如我的這個tif用的是UInt16，也就是無符號的短整型，其他的用的可能就是byte型或其他型別的了。獲取波段資料的資料型別資訊可以如下進行：

cout<<"Data Type:"<<poBand->GetRasterDataType()<<endl;

這裡返回的資訊是一個資料型別的編號，對應列舉型別GDALDataType中的一個值，在此tif的數值是2，對應GDT_UInt16，如果希望返回可讀性比較強的資訊，可以如下

cout<<"Data Type:"<<GDALGetDataTypeName(poBand->GetRasterDataType())<<endl;

返回的是資料型別名的字串：Data Type:UInt16

 

關於波段色彩型別的情況類似資料型別，也是可以由GetColorInterpretation獲取關於型別分類的GDALColorInterp列舉值，同時由GDALGetColorInterpretationName來返回名稱字串：

cout<<"Color interpretation :"<<poBand->GetColorInterpretation()<<endl;

cout<<"Color interpretation:"<<GDALGetColorInterpretationName(poBand->GetColorInterpretation())<<endl;

輸出：

Color interpretation :1

Color interpretation:Gray

表明這個波段是個灰度圖波段。如果是真彩色圖可能得到的結果就是Red，Blue或Green了。

 

GetMaximum和GetMinimum分別返回的是波段資料值可能的最大值和最小值，比如UInt16型別的波段最大值是65535，最小值就是0，知道這兩個值就可以將影象轉換成0～255的數值用於影象的顯示操作了。如：

cout<<"Max :"<<poBand->GetMaximum()<<endl;

cout<<"Min :"<<poBand->GetMinimum()<<endl;

 

終於要說到具體資料的讀寫了。關於波段資料的續寫核心函式就是RasterIO。這個函式可以將影象的某一個子塊讀入或寫入。

CPLErr GDALRasterBand::RasterIO ( GDALRWFlag  eRWFlag, int  nXOff, int  nYOff, int  nXSize, int nYSize, void *  pData, int  nBufXSize, int  nBufYSize, GDALDataType  eBufType, int  nPixelSpace, int  nLineSpace  )

eRWFlag：讀寫標誌位，GF_Read讀，GF_Write寫。

nXOff，nYOff：讀寫資料塊的左上角座標

nXSize，nYSize：讀寫資料塊的xy方向畫素尺寸

pData：讀寫資料塊緩衝區指標

nBufXSize，nBufYSize：緩衝區的畫素寬度和高度

eBufType：資料型別，如GDT_UInt16等

nPixelSpace：pData中一個畫素的位元組大小，值為0的時候就預設是sizeof(eBufType)

nLineSpace ：pData中一行資料的位元組大小，值為0的時候就預設是sizeof(eBufType)*nBufXSize

       比如說我們要讀取第一行資料：

short *pafScanline;

     int   nXSize = poBand->GetXSize();

     pafScanline = ( short *) CPLMalloc(sizeof( short)*nXSize); //分配緩衝區空間

     poBand->RasterIO( GF_Read, 0, 0, nXSize, 1,

     pafScanline, nXSize, 1, GDT_UInt16, 0, 0 );

     CPLFree(pafScanline); //釋放緩衝區空間

 

獲取波段的畫素資料後就可以使用GDI或GL的畫素操作函式繪製影象了。另外，如果希望一次性讀出整個波段的所有資料，根據GDAL文件的說法使用ReadBlock函式效率會更好一些。