Contents
GDAL库的一些细节
关于ReadRaster
缩放
下面是我的一个小例子：我使用了ReadAsArray，返回的是Numeric模块定义的Array，我喜欢它的是因为它排列很工整，看起来很好看。 (这里的数据是用GRASS从spearfish示例数据转出的aspect数据集)
1 >>> import gdal
2 >>> from gdalconst import *
3 >>> dataset = gdal.Open("f:/pyproj/gtif/aspect.tif")
4 >>> band = dataset.GetRasterBand(1)
5 >>> band.ReadAsArray(100,100,5,5,10,10)
6 array([[ 61, 61, 58, 58, 90, 90, 87, 87, 45, 45],
7 [ 61, 61, 58, 58, 90, 90, 87, 87, 45, 45],
8 [ 36, 36, 59, 59, 113, 113, 88, 88, 37, 37],
9 [ 36, 36, 59, 59, 113, 113, 88, 88, 37, 37],
10 [255, 255, 82, 82, 135, 135, 72, 72, 22, 22],
11 [255, 255, 82, 82, 135, 135, 72, 72, 22, 22],
12 [ 45, 45, 129, 129, 144, 144, 255, 255, 255, 255],
13 [ 45, 45, 129, 129, 144, 144, 255, 255, 255, 255],
14 [ 90, 90, 110, 110, 98, 98, 35, 35, 45, 45],
15 [ 90, 90, 110, 110, 98, 98, 35, 35, 45, 45]],'b')
16 >>> band.ReadAsArray(100,100,10,10,10,10)
17 array([[ 61, 58, 90, 87, 45, 255, 255, 117, 65, 50],
18 [ 36, 59, 113, 88, 37, 26, 63, 165, 23, 74],
19 [255, 82, 135, 72, 22, 29, 67, 118, 77, 120],
20 [ 45, 129, 144, 255, 255, 36, 94, 108, 88, 97],
21 [ 90, 110, 98, 35, 45, 88, 109, 121, 92, 73],
22 [ 94, 108, 85, 45, 55, 97, 144, 167, 135, 21],
23 [ 96, 106, 82, 45, 45, 45, 255, 230, 251, 255],
24 [246, 236, 255, 255, 3, 255, 255, 247, 254, 255],
25 [236, 249, 255, 255, 255, 255, 255, 255, 255, 255],
26 [194, 212, 255, 255, 255, 255, 255, 255, 255, 255]],'b')
27 >>> band.ReadAsArray(100,100,20,20,10,10)
28 array([[ 54, 95, 91, 150, 53, 135, 164, 139, 35, 37],
29 [128, 152, 86, 97, 96, 91, 116, 199, 255, 200],
30 [101, 66, 71, 135, 80, 152, 255, 255, 255, 210],
31 [171, 159, 87, 247, 254, 202, 104, 255, 255, 160],
32 [223, 255, 255, 255, 255, 206, 147, 193, 218, 121],
33 [227, 255, 255, 116, 150, 238, 255, 255, 137, 162],
34 [230, 255, 231, 247, 145, 156, 134, 30, 130, 136],
35 [155, 196, 252, 230, 187, 19, 134, 195, 126, 144],
36 [ 80, 54, 85, 105, 163, 140, 192, 95, 154, 146],
37 [150, 83, 71, 161, 173, 255, 255, 120, 149, 180]],'b')
38 >>>
ReadAsArray的原型
看看函数的几个参数的意义。头两个100是取值窗口的左上角在实际数据中所处象元的xy位置。后两个是取值窗口覆盖的区域大小，再后面 是取值窗口取出数组进行缩放后数组的大小。这里需要注意的是这里的buffer大小是根据参数自动分配的，可以不指定，如果不指定，则和第3，4个参数一致。经过5，6两个参数的设置，可以进行缩放。比如真的取值窗口大小可以是20*20，而取出后数组就可以人工设置大小。让他称为10*10的数组，或者是40*40的数组。如果设置20*40则取出的数组对于真实图像来说有了变形。
1 >>> band.ReadAsArray(100,100,10,10)
2 array([[ 61, 58, 90, 87, 45, 255, 255, 117, 65, 50],
3 [ 36, 59, 113, 88, 37, 26, 63, 165, 23, 74],
4 [255, 82, 135, 72, 22, 29, 67, 118, 77, 120],
5 [ 45, 129, 144, 255, 255, 36, 94, 108, 88, 97],
6 [ 90, 110, 98, 35, 45, 88, 109, 121, 92, 73],
7 [ 94, 108, 85, 45, 55, 97, 144, 167, 135, 21],
8 [ 96, 106, 82, 45, 45, 45, 255, 230, 251, 255],
9 [246, 236, 255, 255, 3, 255, 255, 247, 254, 255],
10 [236, 249, 255, 255, 255, 255, 255, 255, 255, 255],
11 [194, 212, 255, 255, 255, 255, 255, 255, 255, 255]],'b')
通 过几个例子可以看到，读取的4个size参数的作用就是把硬盘上指定区域（前四个参数定义）的数据按比例缩放到用户指定区域（后两个定义）内，必要的时候 进行缩放。这里需要注意的是缩的情况（缩的时候是取几个周围点的平均值）如果是调色板数据就可能引发问题（是否可以设置重采样方式我还要再研究一下）。
补：在这里发现了一个帖子，发现RasterIO用的都是最临近(都是最临近？)，而要设置重采样方式只能在BuildPyramid的时候设置了。现在看来BuildPyramid还是有些用的。
范围
现在还有一个问题。就是取值窗口超过实际数据最大的范围怎么办？
1 >>> band.XSize
2 634
3 >>> band.YSize
4 478
5 >>> band.ReadAsArray(630,100,10,10)
6 ERROR 5: Access window out of range in RasterIO(). Requested
7 (630,100) of size 10x10 on raster of 634x478.
8 Traceback (most recent call last):
9 File "<stdin>", line 1, in ?
10 File "D:\Python24\lib\gdal.py", line 837, in ReadAsArray
11 buf_xsize, buf_ysize, buf_obj )
12 File "D:\Python24\lib\gdalnumeric.py", line 178, in BandReadAsArray
13 buf_xsize, buf_ysize, datatype, buf_obj )
14 TypeError: Access window out of range in RasterIO(). Requested
15 (630,100) of size 10x10 on raster of 634x478.
16 >>>
可以看到，出错了。获取数据的时候不能越界，很不好。还要调用的时候去判断越界没。还好用python的Numeric模块可以很方便地扩展矩阵。
效率
对于gtif来说，从横向读取和重纵向读取的效率会相差十分之大（那不是一点的大）！ 可以写一个python文件来验证一下
1 # gtif.py
2 import gdal
3 import time
4 dataset = gdal.Open("f:/gisdata/gtif/zz.tif")
5 band = dataset.GetRasterBand(1)
6 width = dataset.RasterXSize
7 height = dataset.RasterYSize
8 bw = 128
9 bh = 128
10 bxsize = width/128
11 bysize = width/128
12 start = time.time()
13 band.ReadAsArray(0,0,width,height)
14 print time.time()-start
15 start2 = time.time()
16 for i in range(bysize):
17 for j in range(bxsize):
18 band.ReadAsArray(bw*j,bh*i,bw,bh)
19 print time.time()-start2
运行一下得到结果
F:\pyproj>gtif.py 2.35899996758 0.766000032425
然后把最后一段的循环的两个for位置调换一下（缩进要注意），变成：
运行结果是：
F:\pyproj>gtif.py 2.48400020599 24.140999794
天！运行时间是第一次的N倍！切注意！切注意！
另外，我们可以看到，如果把图像分块读取，比不分块读取同样大小的图像效率明显要高出许多。
关于ColorMap
ColorMap的颜色定义
ColorMap的定义在下面有详细的解释 :
颜色表： 一个颜色表可能包含一个或者更多的如下面这种C结构描述的颜色元素。
typedef struct
{
/- gray, red, cyan or hue -/
short c1;
/- green, magenta, or lightness -/
short c2;
/- blue, yellow, or saturation -/
short c3;
/- alpha or blackband -/
short c4;
} GDALColorEntry;
颜色表也有一个色彩解析值。（GDALPaletteInterp）这个值有可能是下面值的一种。并且描述了C1,c2,c3,c4该如何解释。
GPI_GRAY: c1表示灰度值 GPI_RGB: c1表示红，c2表示绿，c3表示蓝，c4表示Alpha GPI_CYMP: c1表示青，c2表示品，c3表示黄，c4表示黑 GPI_HLS: c1表示色调，c2表示亮度，c3表示饱和度
虽然有颜色表示数的区别，但是用GetColorEntry读出的都是4个值，根据PaletteInterp枚举值看截取其中的几个值形成颜色。
ColorMap颜色变动
需要注意的是在gdal使用ColorMap的时候，对原始的ColorMap已经做过变动了。比如geotiff的ColorMap的数据类型是 short，默认的范围是在0～65535，但是在gdal读取出来以后，已经经过了范围压缩。压到了0～255。虽然都是short类型，但是值已经变化。
参考快速开始那张颜色表，可以看到颜色表中的数据是经过（原始数据/65535.0*255）调整过的。 这里在使用的时候可能比较方便，但是如果你是有读取过geotiff原始数据背景的，则需要小心。可能会有习惯性思维的陷阱。因为两者的类型都是short，但是实际的数值不一样（我就曾经以为gdal读取的是geotiff内部的原始ColorMap数据）。
反馈
如果您发现我写的东西中有问题，或者您对我写的东西有意见，请登陆这个论坛。