轉自趙文原文
gdal讀寫圖像分塊處理（精華版）
Review：
用gdal，感覺還不如直接用C++底層函數對遙感數據進行處理。因為gdal進行太多封裝，如果你僅僅只是Geotif等格式進行處理，IO，遍歷，轉換，算法處理等操作，就別用gdal了。如果你想懶省事，那么這篇文章還是或許有些參考價值了。但是不推薦你這么做。
一．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
置１逐像素存儲，置０按波段存儲