使用 RasterIO 读取影像
目录
当我们使用 GDAL 从栅格数据(RasterData)影像中读写数据的时候，我们最常使用的方法就是 RasterIO。 GDAL 在 GDALDataset 类和 GDALRasterBand 类中都提供了 RasterIO 方法。在这两个类中，RasterIO 函数的差别不大。下面我们通过一个实例来具体了解如何使用 RasterIO
原型
在 GDALDataset 中，RasterIO 方法的定义如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
CPLErr GDALDataset::RasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, void * pData, int nBufXSize, int nBufYSize, GDALDataType eBufType, int nBandCount, int * panBandMap, int nPixelSpace, int nLineSpace, int nBandSpace ) |
下面我们来了解各个参数的意义：
eRWFlag
它用来指示 RasterIO 的操作方式，读，还是写。GDALRWFlag 枚举有两个值， GF_Read = 0 为读取数据， GF_Write = 1 为写数据。
nXOff,nYOff
它用来指示 RasterIO 在影像上操作(读或写)范围的起始坐标。在影像中，左上角的坐标为 (0,0)。
nXSize,nYSize
它用来指示 RasterIO 在影像上操作范围的大小。nXSize为宽，nYSize为高。 它以 (nXOff,nYOff)为起点，决定了操作影像的范围的大小。
pData
这是一个指针，指向缓存数据块。如果 eRWFlag 为 GF_Write,函数会将 pData 中的数据写入影像，如果为 GF_Read, 那么函数会将从影像中取得的数据写入到 pData 中。 pData 为 void* 类型，它真实的数据类型将由参数 eBufType 决定。
由于 pData 为 void* 类型，所以它的大小一般使用 字节(byte) 来作单位。它的大小需要最底满足缓冲区的大小。即 缓冲区内存块大小 = 缓冲区长 X 缓冲区宽 X 波段数 X 每波段数据元大小 。用参数来表示：
|
1 2 |
//GDALGetDataTypeSize返回数据类型所占位数，故需要除以 8 int nDataSize = nBufXSize * nBufYSize * nBandCount * GDALGetDataTypeSize(eBufType) / 8 |
nBufXSize,nBufYSize
它用来指定缓冲区的宽高。这两个参数配合 nXSize,nYSize 来实现对影像的缩放。例如，如果 nBufXSize > nXSize,即缓存区宽度比影像操作区宽度大，那么读取数据时，将得到一个在 X 方向放大的缓冲数据。反之则得到缩小的缓冲数据。
值得注意的是，在读取并缩小影像时，GDAL 将自动分析金字塔(Overviews)中数据，选取合适的金字塔层来来读取数据。所以建立合适的金字塔，对读取影像的效率的很大帮助。
eBufType
它也是一个枚举，用来指示操作的数据类型。如果它与 GDALRasterBand 的数据类型不一致，则 RasterIO 将会做类型转换。需要注意的是，将数据从较大的单位向较小的单位进行转换时，RasterIO 会将数据截取而不是比例缩小。例如从 GDT_Int16(Thirty two bit unsigned integer)转为 GDT_Byte(Eight bit unsigned integer)时，超出 255 的部分将会被直接丢弃。
nBandCount,panBandMap
nBandCount为要操作的波段数。panBandMap则为读写的波段的顺序,即我们可以先读取哪一个波段，后读取哪一个波段。它是一个 int 型数组，内容为波段编号，大小为 nBandCount。这个参数可以为 NULL,这样将默认顺序使用影像的前 nBandCount 个波段来读写。需要注意的是，波段编号是从 1 开始的。
nPixelSpace
它用来指示 pData 中的一行(scanline)内，一个像素的起始位置到另一个像素的起始位置间隔的字节数(byte).缺省值 0 表示像素间隔为 eBufType 所占的字节数。
nLineSpace
它用来指示 pData 中每一行(scanline)的起始位置到下一行的起始位置间的字节数。缺省值 0 表示 eBufType 与 nBufXSize 的乘积。
nBandSpace
它用来指示 pData 中，每个波段的起始位置和下一个波段起始位置间的字节数。缺省值 0 表示 nLineSpace 与 nBufYSize 的乘积
实例
下面通过一个实例来深入了解RasterIO的使用。
如图，这是一张 3 波段的 Landsat 影像。宽 6061，高 3907。我们将从影像中读取范围为[(2000,2000),(4000,3500)] 矩形范围内的数据，写入到 1024 X 768 像素的图片中。
在这个例子中，我们使用 cairo 将最终的图像输出到 24bit RGB 格式的 png 图片中。
首先，我们确定了参数 nXOff,nYOff,nXSize,nYSize,nBufXSize,nBufYSize. 我们让 RasterIO 将 2000 * 2500 的影像缩小到 1024 * 768.
然后我们要了解，我们将要写入的图片格式以什么样的方式来存储波段。
在cairo的文档中，我们找到了 CAIRO_FORMAT_RGB24 的说明：
|
1 2 3 |
//CAIRO_FORMAT_RGB24: each pixel is a 32-bit quantity, with // the upper 8 bits unused. Red, Green, and Blue are stored // in the remaining 24 bits in that order. |
可以看出，24位 RGB 格式实际上每个像素占 32 位，高 8 位是无用的，然后由高到低依次为 R,G,B 排列。如下图所示：
如果令像元为 M,像素为 P ,图片宽为 W,高为 H,缓冲区为 BUF ,则有：
P = {MB,MG,MR,MX} ;
BUF = {P0,P1,P2,…,Pn},n=W*H ;
SIZE(BUF) = W * 32 / 8 * H (byte)
我们需要做为，就是控制 panBandMap,nPixelSpace,nLineSpace,nBandSpace 等参数，让 RasterIO 将三个波段的数据按上述顺序写入到缓冲区 pData 中。
- 在该影像中，波段 Band1,Band2,Band3 分别做为 红,绿,蓝 通道，则: panBandMap 为 {3,2,1} 。
- 由前述得到每个像素为 32 位，即: nPixelSpace 为 4
- 每一行大小就是图像的宽，即 nLineSpace 为 4 * 1024
- 波段与波段间的间隔为一个字节，即 nBandSpace 为 1
至此，我们的已经完成了此实例的分析。下面给出完整的实现代码：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 |
#include <iostream> #include <gdal/gdal.h> #include <gdal/gdal_priv.h> #include <cairo/cairo.h>; using namespace std; const char* pszfileName = "E:\\data\\LANDSAT_1_X05.tif"; const char* pszOutfile = "E:\\data\\out.png"; int main(void) { GDALAllRegister(); //init and register all CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //set file name utf8 //open the file to dataset GDALDataset* hDs = (GDALDataset*)GDALOpen(pszfileName, GA_ReadOnly); if (hDs == NULL) return 0; //let data type to GDT_Byte(0~255) GDALDataType dy = GDT_Byte; //init buffer int nBufW = 1024,nBufH = 768; int nSz = nBufW* nBufH * GDALGetDataTypeSize(dy) * 4 / 8; unsigned char* pszData = (unsigned char*)malloc(nSz); memset(pszData, 0, nSz); int nOx = 2000, nOy = 2000, nX = 4000, nY = 3500; int nBandCount = 3; int bandMap[3] = { 3,2,1 }; //do RasterIO CPLErr er = hDs->RasterIO(GF_Read, nOx, nOy, nX - nOx, nY - nOy, pszData, nBufW, nBufH, dy, nBandCount, bandMap, 4, 4 * nBufW , 1); if (er != CE_None) cout << CPLGetLastErrorMsg() << endl; GDALClose(hDs); //write buffer to png file with cairo cairo_surface_t* sfs = cairo_image_surface_create_for_data(pszData, CAIRO_FORMAT_RGB24, nBufW, nBufH , nBufW * 4); cairo_t* ca = cairo_create(sfs); cairo_status_t stat = cairo_surface_write_to_png(sfs, pszOutfile); cairo_surface_destroy(sfs); cairo_destroy(ca); return 0; } |
最终我们得到的图片如下：
扩展
在此例中，因为恰好 unused 在高 8 位，所以我们很方便地将 panBandMap 置为 {3,2,1}. 将 nPixelSpace 置为 4 ，每个像素前 3 个字节先写入三个波段，然后空出 1 个字节来做为 unused 。 但是，如果 unused 不在高 8 位，在高位第 2 个字节呢？ 此时，pData中,每 4 个字节的排列为：
P = {MB,MG,MX,MR}
我们将如何向 RasterIO 传入参数呢？ 解决方案有多种：
方案 1
在分配缓冲区空间时，多申请3个字节的空间。将 pData[ 0 ] 做为起始缓冲区的起始地址传入 RasterIO. panBandMap为 {1,3,2}, 让 RasterIO 按 RBGX 的顺序写入缓冲区。然后将最后3个字节的内容(BGX)移到pData[ 0 ],此时的 pData 内每个像素的排列则为 BGXR。
方案 2
令 nBandCount 为 4 ，panBandMap 为 {3,2,1,1}, 其中，panBandMap[ 0 ]为占位波段，在RasterIO 完成后，将每个像素内此波段的置 0
方案 3
前面我们提到，不仅 GDALDataset 为提供了 RaseterIO，GDALRasterBand 也提供了 RasterIO, 用于单独读写每个波段的数据。我们可以将每个波段分别写入不同的缓冲区，然后就可以按任意方式进行组合
这种方案比前几种稍显复杂，但更具有灵活性，是一种常用的方式。