<> = 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') >>> }}} !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 "", 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|这个论坛]]。