5.2. GDAL库的一些细节#
== 关于ReadRaster ==
=== 缩放 ===
下面是我的一个小例子：我使用了!ReadAsArray，返回的是Numeric模块定义的Array，我喜欢它的是因为它排列很工整，看起来很好看。 (这里的数据是用GRASS从spearfish示例数据转出的aspect数据集)
#! python
>>> import gdal
>>> from gdalconst import *
>>> dataset = gdal.Open("f:/pyproj/gtif/aspect.tif")
>>> band = dataset.GetRasterBand(1)
>>> band.ReadAsArray(100,100,5,5,10,10)
array([[ 61, 61, 58, 58, 90, 90, 87, 87, 45, 45],
[ 61, 61, 58, 58, 90, 90, 87, 87, 45, 45],
[ 36, 36, 59, 59, 113, 113, 88, 88, 37, 37],
[ 36, 36, 59, 59, 113, 113, 88, 88, 37, 37],
[255, 255, 82, 82, 135, 135, 72, 72, 22, 22],
[255, 255, 82, 82, 135, 135, 72, 72, 22, 22],
[ 45, 45, 129, 129, 144, 144, 255, 255, 255, 255],
[ 45, 45, 129, 129, 144, 144, 255, 255, 255, 255],
[ 90, 90, 110, 110, 98, 98, 35, 35, 45, 45],
[ 90, 90, 110, 110, 98, 98, 35, 35, 45, 45]],'b')
>>> band.ReadAsArray(100,100,10,10,10,10)
array([[ 61, 58, 90, 87, 45, 255, 255, 117, 65, 50],
[ 36, 59, 113, 88, 37, 26, 63, 165, 23, 74],
[255, 82, 135, 72, 22, 29, 67, 118, 77, 120],
[ 45, 129, 144, 255, 255, 36, 94, 108, 88, 97],
[ 90, 110, 98, 35, 45, 88, 109, 121, 92, 73],
[ 94, 108, 85, 45, 55, 97, 144, 167, 135, 21],
[ 96, 106, 82, 45, 45, 45, 255, 230, 251, 255],
[246, 236, 255, 255, 3, 255, 255, 247, 254, 255],
[236, 249, 255, 255, 255, 255, 255, 255, 255, 255],
[194, 212, 255, 255, 255, 255, 255, 255, 255, 255]],'b')
>>> band.ReadAsArray(100,100,20,20,10,10)
array([[ 54, 95, 91, 150, 53, 135, 164, 139, 35, 37],
[128, 152, 86, 97, 96, 91, 116, 199, 255, 200],
[101, 66, 71, 135, 80, 152, 255, 255, 255, 210],
[171, 159, 87, 247, 254, 202, 104, 255, 255, 160],
[223, 255, 255, 255, 255, 206, 147, 193, 218, 121],
[227, 255, 255, 116, 150, 238, 255, 255, 137, 162],
[230, 255, 231, 247, 145, 156, 134, 30, 130, 136],
[155, 196, 252, 230, 187, 19, 134, 195, 126, 144],
[ 80, 54, 85, 105, 163, 140, 192, 95, 154, 146],
[150, 83, 71, 161, 173, 255, 255, 120, 149, 180]],'b')
>>>
}}}
5.2.1. ReadAsArray的原型#
- ::
#! python >>> help(band.ReadAsArray) Help on method ReadAsArray in module gdal: ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None , buf_ysize=None, buf_obj=None) method of gdal.Band instance
看看函数的几个参数的意义。头两个100是取值窗口的左上角在实际数据中所处象元的xy位置。后两个是取值窗口覆盖的区域大小，再后面 是取值窗口取出数组进行缩放后数组的大小。这里需要注意的是这里的buffer大小是根据参数自动分配的，可以不指定，如果不指定，则和第3，4个参数一致。经过5，6两个参数的设置，可以进行缩放。比如真的取值窗口大小可以是20*20，而取出后数组就可以人工设置大小。让他称为10*10的数组，或者是40*40的数组。如果设置20*40则取出的数组对于真实图像来说有了变形。
#! python
>>> band.ReadAsArray(100,100,10,10)
array([[ 61, 58, 90, 87, 45, 255, 255, 117, 65, 50],
[ 36, 59, 113, 88, 37, 26, 63, 165, 23, 74],
[255, 82, 135, 72, 22, 29, 67, 118, 77, 120],
[ 45, 129, 144, 255, 255, 36, 94, 108, 88, 97],
[ 90, 110, 98, 35, 45, 88, 109, 121, 92, 73],
[ 94, 108, 85, 45, 55, 97, 144, 167, 135, 21],
[ 96, 106, 82, 45, 45, 45, 255, 230, 251, 255],
[246, 236, 255, 255, 3, 255, 255, 247, 254, 255],
[236, 249, 255, 255, 255, 255, 255, 255, 255, 255],
[194, 212, 255, 255, 255, 255, 255, 255, 255, 255]],'b')
通 过几个例子可以看到，读取的4个size参数的作用就是把硬盘上指定区域（前四个参数定义）的数据按比例缩放到用户指定区域（后两个定义）内，必要的时候 进行缩放。这里需要注意的是缩的情况（缩的时候是取几个周围点的平均值）如果是调色板数据就可能引发问题（是否可以设置重采样方式我还要再研究一下）。
补：在[[http://lists.maptools.org/pipermail/gdal-dev/2005-May/005646.html|这里]]发现了一个帖子，发现RasterIO用的都是最临近(都是最临近？)，而要设置重采样方式只能在!BuildPyramid的时候设置了。现在看来!BuildPyramid还是有些用的。
=== 范围 ===
现在还有一个问题。就是取值窗口超过实际数据最大的范围怎么办？
#! python
>>> band.XSize
634
>>> band.YSize
478
>>> band.ReadAsArray(630,100,10,10)
ERROR 5: Access window out of range in RasterIO(). Requested
(630,100) of size 10x10 on raster of 634x478.
Traceback (most recent call last):
File "<stdin>", line 1, in ?
File "D:\Python24\lib\gdal.py", line 837, in ReadAsArray
buf_xsize, buf_ysize, buf_obj )
File "D:\Python24\lib\gdalnumeric.py", line 178, in BandReadAsArray
buf_xsize, buf_ysize, datatype, buf_obj )
TypeError: Access window out of range in RasterIO(). Requested
(630,100) of size 10x10 on raster of 634x478.
>>>
可以看到，出错了。获取数据的时候不能越界，很不好。还要调用的时候去判断越界没。还好用python的Numeric模块可以很方便地扩展矩阵。
=== 效率 ===
对于gtif来说，从横向读取和重纵向读取的效率会相差十分之大（那不是一点的大）！ 可以写一个python文件来验证一下
#! python
# gtif.py
import gdal
import time
dataset = gdal.Open("f:/gisdata/gtif/zz.tif")
band = dataset.GetRasterBand(1)
width = dataset.RasterXSize
height = dataset.RasterYSize
bw = 128
bh = 128
bxsize = width/128
bysize = width/128
start = time.time()
band.ReadAsArray(0,0,width,height)
print time.time()-start
start2 = time.time()
for i in range(bysize):
for j in range(bxsize):
band.ReadAsArray(bw*j,bh*i,bw,bh)
print time.time()-start2
运行一下得到结果
F:\pyproj>gtif.py
2.35899996758
0.766000032425
}}}
然后把最后一段的循环的两个for位置调换一下（缩进要注意），变成：
- ::
#! python for j in range(bxsize):
- for i in range(bysize):
band.ReadAsArray(bw*j,bh*i,bw,bh)
运行结果是：
- ::
F:pyproj>gtif.py 2.48400020599 24.140999794
天！运行时间是第一次的N倍！切注意！切注意！
另外，我们可以看到，如果把图像分块读取，比不分块读取同样大小的图像效率明显要高出许多。
== 关于ColorMap ==
=== ColorMap的颜色定义 ===
!ColorMap的定义在[[http://www.gdal.org/gdal_datamodel.html|下面]]有详细的解释 :
颜色表： 一个颜色表可能包含一个或者更多的如下面这种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类型，但是值已经变化。
参考[[http://wiki.woodpecker.org.cn/moin/lilin/gdal-introduce#head-ddf209b471ede52abec271c974ff6e89c5daed49|快速开始]]那张颜色表，可以看到颜色表中的数据是经过（原始数据/65535.0*255）调整过的。 这里在使用的时候可能比较方便，但是如果你是有读取过geotiff原始数据背景的，则需要小心。可能会有习惯性思维的陷阱。因为两者的类型都是short，但是实际的数值不一样（我就曾经以为gdal读取的是geotiff内部的原始!ColorMap数据）。
=== 反馈 ===
如果您发现我写的东西中有问题，或者您对我写的东西有意见，请登陆[[http://groups.google.com/group/geosings?hl=zh-CN|这个论坛]]。