/* 版權所有者：趙文   */

一．gdal進行資料操作
在安裝好gdal後，即可呼叫gdal庫中的函式。
（需要包含的標頭檔案：gdal_priv.h）
1．開啟資料集
使用gdal庫進行資料（影像）操作的第一步就是開啟一個數據集。對於“資料集”這個名詞大家可能不會太習慣，但是對於一般的格式來說，一個“資料集”就是一個檔案，比如一個TIFF檔案就是一個以tiff為副檔名的檔案。但是對於眾多RS資料來說，一個數據集包含的絕對不僅僅是一個檔案。對於很多RS資料，他們把一張影象分成數個影象檔案，然後放在一個資料夾中，用一些額外的檔案來組織它們之間的關係，形成一個“資料集”(有點難以理解，暫且放過)。下面我們由給定的檔案路徑檔名開啟一個tiff(GeoTIFF)檔案。（任何支援的格式，開啟方式都是這樣）

CString strFilePath;
StrFilePath=’d:/rsdata/2005_234.tif’;
GDALDataSet *poDataset; //GDAL資料集
GDALAllRegister();
poDataset = (GDALDataset *) GDALOpen(strFilePath, GA_ReadOnly );

這樣我們就打開了這個檔案。通過資料集poDataset即可呼叫各功能函式，如:
GetRasterCount();//獲取影象波段數；
GetRasterXSize();//獲取影象寬度
GetRasterYSize();//獲取影象高度
GetRasterBand();//獲取影象某一波段
GetGeoTransform(double *p);//獲取影象地理座標資訊長度為六的陣列
RasterIO();//對影象資料進行縮放讀和寫
……
(更具體的API列表可以看這裡)。

2．獲取影象資訊、讀取影象
開啟檔案後，下面要做的就是獲取檔案的相關資訊儲存在相應變數中，將影象資料讀入記憶體中，等待後續處理了。

2．1 獲取基本資訊
因為不同格式資料所包含的相關資訊有所不同，一般情況下我們需要得到影象的高、寬、波段數、地理座標資訊，資料型別等。Gdal庫中有相應的函式實現這些功能。下面的程式碼實現獲取這些資訊：

int nBandCount=poDataset->GetRasterCount();
int   nImgSizeX=poDataset->GetRasterXSize();
int     nImgSizeY=poDataset->GetRasterYSize();
double        adfGeoTransform[6];
poDataset->GetGeoTransform( adfGeoTransform );
//如果影象不含地理座標資訊，預設返回值是：（0、1、0、0、0、1），表中第四列表示了帶//有地理座標資訊的資料格式
// adfGeoTransform[0]是左上角像元的東座標；
// adfGeoTransform[3]是左上角像元的北座標；
// adfGeoTransform[1]是像元寬度；
// adfGeoTransform[5]是像元高度；
//影象其他點座標計算公式：
//Xp = adfGeoTransform [0] + P*adfGeoTransform [1]+L*adfGeoTransform [2];
//Yp = adfGeoTransform [3] + P*adfGeoTransform [4] + L*adfGeoTransform [5];
//注：用GetGeoTransform()並不能十分合理的表示影象地理座標，當影像範圍很大
//時，這種座標表示方法將不適用。

2．2 將影象資料按照要求讀入記憶體
影象的讀寫是通過RasterIO()實現的。RasterIO()功能十分強大,它可以把影象上指定大小的矩形象素塊以縮放的形式按指定的資料型別輸出或輸入到使用者指定大小的緩衝區中。
原型：
CPLErr GDALDataset::RasterIO(
GDALRWFlag eRWFlag,    //讀寫標記如果為GF_Read，則是將影像內容寫入記憶體，如果
//為GF_Write，則是將記憶體中內容寫入檔案。
int nXOff, int nYOff, //相對於影象左上角頂點(從零開始)的行列偏移量
int nXSize, int nYSize,    //要讀寫的塊在x方向的象素個數和y方向的象素列數
void * pData,                //指向目標緩衝區的指標，由使用者分配
int nBufXSize, int nBufYSize,//目標塊在x方向上和y方向上的大小
GDALDataType eBufType,   //目標緩衝區的資料型別，原型別會自動轉換為目標型別
int nBandCount,           //要處理的波段數
int * panBandMap,        //記錄要操作的波段的索引(波段索引從1開始)的陣列，若為空
//則陣列中存放的是前nBandCount個波段的索引
int nPixelSpace,        //X方向上兩個相鄰象素之間的位元組偏移，預設為0,
//則列間的實際位元組偏移由目標資料型別eBufType確定
int nLineSpace,     //y方向上相鄰兩行之間的位元組偏移, 預設為0,則行間的實際位元組
//偏移為eBufType * nBufXSize
int nBandSpace //相鄰兩波段之間的位元組偏移，預設為0,則意味著波段是順序結構
//的，其間位元組偏移為nLineSpace * nBufYSize
);
下面將通過例項演示其使用方法,實現的是將7波段影象中的第2 3 4波段按照3 4 2的順序讀入記憶體中,逐畫素儲存:
   int bandmap[7];
   bandmap[0]=3;
   bandmap[1]=4;
   bandmap[2]=2;
   Scanline=(nImagSizex*8+31)/32*4;
　BYTE   pafScan=( BYTE )CPLMalloc(sizeof(byte) *nImgSizeX*nImgSizeY*３)
　poDataset->RasterIO( GF_Read, 0, 0,nImgSizeX,nImgSizeY, pafScan,
nImgSizeX,nImgSizeY,GDT_Byte,３,bandmap,３, Scanline*３,1 );

以這種方式讀取之後，直接可構建點陣圖進行顯示。這裡可以按照自己的需要進行其他方式讀取。以上讀取方式僅僅為了顯示方便，如進行影象處理相關運算，則按波段全部讀出會比較方便:
poDataset->RasterIO( GF_Read, 0, 0,nImgSizeX,nImgSizeY, pafScan,
nImgSizeX,nImgSizeY,GDT_Byte,bandcount,0,0,0,0);
之前開闢記憶體改為：
BYTE   pafScan=new byte[nImgSizeX*nImgSizeY*bandcount];

將影象資料讀入記憶體後，即可通過指標pafScan對影象進行你想要進行的操作了。

3．另存影象
另存檔案其實就是先建立一個新的檔案，然後將資料寫入新檔案中。
使用gdal建立新檔案有兩種方式：Create()和CreateCopy();有些檔案格式支援Create()函式（見第一頁表格第三列），可以使用Create（）直接建立此類格式的檔案，而其他不支援Create（）函式的影象格式，需要先建立tiff格式檔案，然後複製生成目標格式檔案。
CString strFilePath1;//輸入影象檔案路徑名
CString strFilePath2;//另存為的tiff格式影象路徑
CString strFilePath2; //另存為的jpeg格式影象路徑
GDALDataset *poDataset1; //GDAL資料集
GDALDataset *poDataset2; //待建立的GDAL資料集
GDALDataset *poDataset2; //待建立的GDAL資料集
GDALDriver *poDriver;    //驅動，用於建立新的檔案
GDALAllRegister();
poDataset1 = (GDALDataset *) GDALOpen(strFilePath1, GA_ReadOnly );
// 開啟檔案
if( poDataset1 == NULL )
{
AfxMessageBox("檔案開啟失敗!!!");
m_FileFlag=FALSE;
return;
}
//Create方式建立新的tiff檔案：
nBandCount=poDataset1->GetRasterCount();
nImgSizeX=poDataset1->GetRasterXSize();
nImgSizeY=poDataset1->GetRasterYSize();
//獲取影像相關資訊

//建立新檔案
Cstring fomat；
fomat="GTiff"
poDriver = GetGDALDriverManager()->GetDriverByName(fomat);
//設定檔案型別，表格第二列為格式標識，儲存為其他格通過改變fomat值即可

//獲取格式型別
char **papszMetadata = poDriver->GetMetadata();

//根據檔案路徑檔名，影象寬，高，波段數，資料型別，檔案型別，建立新的資料集
poDataset2=poDriver->Create(strFilePath2,nImgSizeX,nImgSizeY,nBandCount,GDT_Byte,papszMetadata);
//座標賦值
double        adfGeoTransform[6]={0，1，0，0，0，1};
poDataset2->SetGeoTransform(adfGeoTransform);

//將原影象資料讀出，進行相應處理後，寫入新檔案
BYTE *ppafScan= new BYTE [nImgSizeX * nImgSizeY *nBandCount];
poDataset1->RasterIO( GF_Read,0,0, nImgSizeX, nImgSizeY, ppafScan,
nImgSizeX, nImgSizeY, GDT_Byte,nBandCount,0,0, 0,0 );
      .
      .
『中間對影象進行處理運算』
         .
//將影象資料寫入新影象中
poDataset2->RasterIO(GF_Write, 0,0,   nImgSizeX, nImgSizeY,
ppafScan,pafsizex,pafsizey, GDT_Byte,nBandCount,0,0, 0,0 );

delete poDataset2;
delete poDriver;
//影象另存完畢

CreateCopy方式建立jpeg格式檔案：
接上面的過程，先不delete，（即已經完成用create方式先將運算完畢的影象建立為tiff格式）
fomat="Jpeg"
poDriver = GetGDALDriverManager()->GetDriverByName(fomat);
poDataset3=poDriver->CreateCopy(strFilePath3,poDataset2,1,papszMetadata,NULL,NULL);
delete poDataset3;
delete poDataset2;
delete podriver;


二．        使用RasterIO()對大影象進行分塊操作

RasterIO（）函式能夠對影象任意指定區域任意波段的資料按指定資料型別，指定排列方式讀入記憶體和寫入檔案中，因此可以實現對大影像的分塊讀，運算，寫操作。對於大影象處理，按照傳統方法，首先要將影象所有資料讀入記憶體中，進行相應操作後，再一次性將處理好的資料寫入檔案中，這樣需要耗費很大記憶體，容易記憶體溢位，而且存續可執行行差。採用分塊處理技術，一幅1G的影像，在整個資料處理過程中，可以只佔用幾十兆的記憶體，而且運算量不會增加。下面通過一個示例加以演示：
/所有波段分塊處理示例
void CTestzwDoc::OnLowers()
{
Inoutput dlg; //獲取檔案路徑的對話方塊類
if (dlg.DoModal()==IDCANCEL)
{
return;
}
CString strFilePath1(dlg.m_input);
CString strFilePath2(dlg.m_output);
GDALDataset *poDataset1; //GDAL資料集
GDALDataset *poDataset2; //GDAL資料集
GDALDriver *poDriver;
GDALAllRegister();

poDataset1 = (GDALDataset *) GDALOpen(strFilePath1, GA_ReadOnly );
if( poDataset1 == NULL )
{
AfxMessageBox("檔案開啟失敗!!!");
m_FileFlag=FALSE;
return;
}
int BandCount=poDataset1->GetRasterCount();
int nImgSizeX=poDataset1->GetRasterXSize();
int nImgSizeY=poDataset1->GetRasterYSize();

//建立新檔案
CString format;
format="Gtiff";
poDriver = GetGDALDriverManager()->GetDriverByName(format);
char ** papszMetadata = poDriver->GetMetadata();
poDataset2=poDriver->Create(strFilePath2,nImgSizeX,nImgSizeY,nBandCount,GDT_Byte,papszMetadata);
//設定影象座標，缺少這一步，建立的影象用erdas開啟將無法正常現實
poDataset2->SetGeoTransform(adfGeoTransform );
/////////////////////////////////////////////////////
//分塊處理.將影像分成很多512*512大小的塊，通過迴圈對每一塊進行處理
int nxNum=(nImgSizeX-1)/512+1;//計算列方向上塊數
int nyNum=(nImgSizeY-1)/512+1;//計算行方向塊數
int pafsizex;                    //當前塊寬度
int pafsizey;                    //當前塊高度
BYTE * lp;
BYTE *ppafScan= new BYTE [512*512*nBandCount];
for (int nYI=0;nYI<nyNum;nYI++)
for (int nXI=0;nXI<nxNum;nXI++)
{  
      pafsizex=512;
       pafsizey=512;
       //行列末尾小塊處理
       if (nXI==nxNum-1)pafsizex=(nImgSizeX-1)%512+1;
       if (nYI==nyNum-1)pafsizey=(nImgSizeY-1)%512+1;
      //讀取當前塊資料
       poDataset1->RasterIO( GF_Read, nXI*512, nYI*512,pafsizex,
pafsizey,ppafScan,pafsizex,pafsizey,GDT_Byte,BandCount,0,0,0,0）;
       //對當前塊進行處理，邊緣提取
       for(int nnum=0;nnum<nBandCount;nnum++)
       for (int i=0;i<pafsizey;i++)
          for(int j=0;j<pafsizex;j++)
       {         
       {
          lp=ppafScan+nnum*pafsizex*pafsizey+i*pafsizex+j;
         if(j==pafsizex-1&&i!=pafsizey-1)
            *lp=abs(*lp-*(lp+pafsizex));
          else if (i==pafsizey-1&&j==pafsizex-1)
            *lp=0;
          else if (i==pafsizey-1&&j!=pafsizex-1)
            *lp=abs(*lp-*(lp+1));
else                                *lp=abs((*lp)-*(lp+1))+abs(*lp-*(lp+pafsizex));
//邊緣處理是難點
         if (*lp<15)
            *lpp=0;
          else if(15<=*lpp<30)
            *lpp=128;
         else if(*lpp>=30)
            *lpp=255;
       }
   }
       //將當前塊資料寫入新影象相應位置
         poDataset2->RasterIO(GF_Write,nXI*512,nYI*512,pafsizex,pafsizey,ppafScan,pafsizex,pafsizey, GDT_Byte,nBandCount,0,0, 0,0 );
}
   delete []ppafScan;
//寫操作完成後必須釋放，不然寫入操作不成功
delete poDataset2;
delete poDataset1;
delete poDriver;
delete dlg;
}

在前面一篇介紹gdal庫讀取和儲存影象的文章中，有很多不足之處，個人覺得其精華在於在記憶體中建立點陣圖，並進行快速顯示部分。
兩個相鄰象素之間的位元組偏移，則如按波段儲存，則置為0
行偏移量，如按波段儲存，則置為0
置１逐畫素儲存，置０按波段儲存