gdal读写图像分块处理（精华版）

Posted 2020-07-01 Libra

一．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
置１逐像素存储，置０按波段存储