使用 RasterIO 读取影像

Posted on Edited on

gdalbaner

当我们使用 GDAL 从栅格数据(RasterData)影像中读写数据的时候，我们最常使用的方法就是 RasterIO。 GDAL 在 GDALDataset 类和 GDALRasterBand 类中都提供了 RasterIO 方法。在这两个类中，RasterIO 函数的差别不大。下面我们通过一个实例来具体了解如何使用 RasterIO

原型[

](https://code.wandoer.com/wp-content/uploads/2015/06/gdalbaner.png)

在 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 像素的图片中。

[

raster

](https://code.wandoer.com/wp-content/uploads/2016/06/raster-1.jpg)

在这个例子中，我们使用 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 排列。如下图所示：

[

rgb_order

](https://code.wandoer.com/wp-content/uploads/2016/06/rgb_order-1.png)

如果令像元为 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;
}

最终我们得到的图片如下：

[

out

](https://code.wandoer.com/wp-content/uploads/2016/06/out.png)

扩展

在此例中，因为恰好 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, 用于单独读写每个波段的数据。我们可以将每个波段分别写入不同的缓冲区，然后就可以按任意方式进行组合
这种方案比前几种稍显复杂，但更具有灵活性，是一种常用的方式。