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的原型

   1 >>> help(band.ReadAsArray)
   2 Help on method ReadAsArray in module gdal:
   3 ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None
   4 , 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则取出的数组对于真实图像来说有了变形。

   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位置调换一下(缩进要注意),变成:

   1 for j in range(bxsize):
   2     for i in range(bysize):
   3         band.ReadAsArray(bw*j,bh*i,bw,bh)

运行结果是:

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数据)。

反馈

如果您发现我写的东西中有问题,或者您对我写的东西有意见,请登陆这个论坛

lilin/gdal-notice (last edited 2009-12-25 07:10:13 by localhost)