您的位置:首页 > 其它

gdal使用经验(一)影像读取-生成金字塔

2014-12-20 15:34 260 查看
1.关于rasterio

################################################################

对于http://www.gdal.org/gdal_tutorial.html的翻译

读取栅格数据。

有几种方法来读取栅格数据。但是最普遍的就是使用GDALRasterBand::RasterIO() 。这个方法自动会照顾到数据类型转换,以及取数据的时候对数据窗口的放大缩小。下面是几个例子可以说明在图像中如何读取第一行数据到一个简单得设定了大小的数据内存缓冲中。其中有的例子有把数据转换成浮点类型。

C++:

        float *pafScanline;

        int   nXSize = poBand->GetXSize();

        pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize);

        poBand->RasterIO( GF_Read, 0, 0, nXSize, 1, 

                          pafScanline, nXSize, 1, GDT_Float32, 

                          0, 0 );

 C:

        float *pafScanline;

        int   nXSize = GDALGetRasterBandXSize( hBand );

        pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize);

        GDALRasterIO( hBand, GF_Read, 0, 0, nXSize, 1, 

                      pafScanline, nXSize, 1, GDT_Float32, 

                      0, 0 );

 Python:

        scanline = band.ReadRaster( 0, 0, band.XSize, 1, /

                                     band.XSize, 1, GDT_Float32 )

注意在python的例子中,返回的数据被组织成字符串形式。而且返回的数据大小是Xsize*4字节的未加工的二进制浮点数据。你可以用python的标准包中的struct模块来把它转化成python的数据。

代码

     import struct

        tuple_of_floats = struct.unpack('f' * b2.XSize, scanline)

这样返回的tuple_of_floats就是一个经过unpack拆分后的float型数组。(译者)

RasterIO的具体形式如下:

CPLErr GDALRasterBand::RasterIO( GDALRWFlag eRWFlag,

                                 int nXOff, int nYOff, int nXSize, int nYSize,

                                 void * pData, int nBufXSize, int nBufYSize,

                                 GDALDataType eBufType,

                                 int nPixelSpace,

                                 int nLineSpace )

注意:这个方法可以同时用于读和写。你可以使用RWFlag来设置(或者GF_Read或者GF_Write)读写权限。nXOff,nYOff,nXSize,nYSize参数用来描述栅格在硬盘上读取或者写入数据的窗口大小。该窗口不一定要落在边界上,虽然这样存取可能更加有效。

pData是指向缓冲数据的指针。如果是写入,pData就是要写到实际数据中的数据,如果是读取,pData就是从实际数据中获取的数据值存放的地方。它的真实的数据类型是什么是通过eBufType参数设定的。如GDT_Float32,或者是GDT_Byte。RasterIo会照顾到从缓冲到实际数据的转换。注意:当从浮点数转换到整形时,数据将四舍五入。而且从一个过大的存储单位转到较小的存储单位时所进行的操作是截断而不是按比例缩小所有的数,比如原来实际数据中float中存取的数据范围超过了缓冲数据指定数据类型byte的最大范围,操作将把超过byte的范围截去顶端数据,而不是把每个float缩小到0~256范围内。

nBufXSize和nBufYSize值描述缓冲区的大小。当加载的数据是完整分辨率的时候(原始数据,没有缩放)它们应该设成和取值窗口一样的大小。但是在加载时使用了缩小分辨率。它们需要设置得比取值窗口小。在这种情况下,RasterIO利用略缩图组(overviews)中的某个合适的略缩图来做存取将更加有效。

nPixelSpace和nLineSpace一般情况下是0作为缺省值。但是,它们可以用于控制用于存取的内存数据缓冲。可以利用它们把包含的象元数据按另一种组织形式读入内存缓冲内。

################################################################
http://www.gdal.org/classGDALRasterBand.html#1eef7a21cf4809425c3edced99aa088f中的解释:
这个方法允许读取GDALRasterBand的一个区域的数据到一个缓冲或者把一个缓冲中的数据写入GDALRasterBand的指定区域内。它自动考虑到了在波段数据和指定缓冲之间如果用eBufType指定数据类型后不同的数据类型间的转换。它也考虑到了如果缓冲大小(nBufXSize * nBufYSize)和实际数据读取窗口大小(nXSize * nYSize)不一样的时候的数据替代。

nPixelSpace and nLineSpace参数允许读取或者写入非常规组织的缓冲数据。这首先可以用于一个缓冲中包含多个波段数据,并且各个数据之间是交叉排列的。(比如一个数据表中的数据组织是RGBRGBRGB……而普通的数据可能是RRR……GGG……BBB……)

一些数据格式如果从低分辨率图中读取数据可能会更有效率。

在全分辨率中存取的最快的操作是使用 GetBlockSize()获得波段范围, 然后用 ReadBlock()和 WriteBlock()存取数据。

这个方法和C中的GDALRasterIO是一样的。

参数:

    eRWFlag    可选GF_Read来读取数据,GT_Write来写入数据

    nXOff    读取区域的左上角象元在整幅图中偏移左边框的象元数

    nYOff    读取区域的左上角象元在整幅图中偏移上边框的象元数

    nXSize    要存取的实际区域的宽度

    nYSize    要存取的实际区域的高度

    pData    存取数据的数据缓冲。大小至少为nBufXSize*nBufYSize*eBufType类型表示的大小。数据组织顺序是从左到右,从上到下。间隔由nPixelSpace和nLineSpace来控制。

    nBufXSize    缓冲区域的宽

    nBufYSize    缓冲区域的高

    eBufType    缓冲中的数据类型。在需要的时候会和实际波段中的数据类型自动进行互相转换。

    nPixelSpace    在一个扫描行中的一个象元的字节偏移起始点到下一个象元字节偏移起始点之间的字节间隔。如果使用默认的0,则使用eBufType作为实际的两个象元其实字节之间的字节间隔。

    nLineSpace    在一行数据和下一行数据之间的起始字节间的间隔。如果使用偏移值(是不是为0?译者)则用eBufType*nBufXSize来表示实际间隔。

返回值:

    CE_Failure是存取失败,否则返回CE_None

#################################################################

下面是我的一个小例子:我使用了python中读取数据的另一种方式,使用的是ReadAsArray,作用等同于RasterIO,返回的是Numeric的Array,这样就不要用struct去拆分二进制字符串了。而且通过Numeric的数组操作,很容易看似很专业地变换整个矩阵(这是一个C的python扩展,非常快!)。但是我喜欢它的是因为它排列很工整,看起来很好看。

>>> import gdal

>>> from gdalconst import *

>>> dataset = gdal.Open("f:/pyproj/gtif/asp.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')

>>>

可以看出利用RasterIO的几个参数的意义。头两个100是取值窗口的左上角在实际数据中所处象元的xy位置。后两个是取值窗口覆盖的区域大小,再后面是buffer的大小。这里需要注意的是这里的buffer大小是根据参数自动分配的,可以不指定,如果不指定,则和第3,4个参数一致。而不是像C/C++那样需要手工分配。

>>> 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')

ReadAsArray的原型

>>> 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

令人惊奇的是……居然可以不要任何参数!不过输出当然也什么都没有!

通过几个例子可以看到,读取的4个size参数的作用就是把硬盘上指定区域(前四个参数定义)的数据按比例缩放到用户指定区域(后两个定义)内,必要的时候进行缩放。这里需要注意的是缩的情况(缩的时候是取几个周围点的平均值)如果是调色板数据就可能引发问题(是否可以设置重采样方式我还要再研究一下)。所以还是用BuildOverviews来进行一下缩放。然后获取吧。

现在还有一个问题。就是取值窗口超过实际数据最大的范围怎么办?

>>> 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可以很方便地扩展矩阵。但是C/C++就不得不用for来解决了(或者有没有更好的解决方法,我再看看)。

如果是读取全分辨率的图像,可以试试ReadBlock,看看ReadBlock在越界时会有什么反映。但是偏偏python没有,-_-!

>>> dir(band)

['AdviseRead', 'Checksum', 'ComputeBandStats', 'ComputeRasterMinMax', 'DataType'

, 'Fill', 'FlushCache', 'GetDescription', 'GetHistogram', 'GetMaximum', 'GetMeta

data', 'GetMinimum', 'GetNoDataValue', 'GetOffset', 'GetOverview', 'GetOverviewC

ount', 'GetRasterColorInterpretation', 'GetRasterColorTable', 'GetScale', 'ReadA

sArray', 'ReadRaster', 'SetDescription', 'SetMetadata', 'SetNoDataValue', 'SetRa

sterColorInterpretation', 'SetRasterColorTable', 'WriteArray', 'WriteRaster', 'X

Size', 'YSize', '__doc__', '__init__', '__module__', '_o']

>>>

没辙!(我可不想写冗长的C++代码去验证)

最后还要尤其注意一点!!!!!!!(尤其应该十分注意)

那就是RasterIO的读取效率

对于gtif来说,从横向读取和重纵向读取的效率会相差十分之大(那不是一点的大)!

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

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倍!切注意!切注意!

 

2、关于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表示饱和度

颜色表的作用就是把一个颜色和一个栅格象元联系起来,象元的值对应的是一个颜色表的下标。这意味着颜色总是以0为最小然后开始递增。没有规定在颜色表中查找对应颜色之前需要指定一个预缩放比例因子的机制(There is no provision for indicating a prescaling mechanism before looking up in the color table这句没看太懂)。

##################################################################

需要注意的是在gdal使用ColorMap的时候,对原始的ColorMap已经做过变动了。比如geotiff的ColorMap的数据类型是 short,默认的范围是在0~65535,但是在gdal读取出来以后,已经经过了范围压缩。压到了0~255。虽然都是short类型,但是值已经变化

>>> colortab = band.GetRasterColorTable()

>>> dir(colortab)

['Clone', 'GetColorEntry', 'GetColorEntryAsRGB', 'GetCount', 'GetPaletteInterpre

tation', 'SetColorEntry', '__del__', '__doc__', '__init__', '__module__', '__str

__', '_o', 'own_o', 'serialize']

>>> ct = []

>>> for i in range(colortab.GetCount()):

...     ct.append(colortab.GetColorEntry(i))

...

>>> ct

[(0, 0, 0, 255), (0, 0, 28, 255), (0, 0, 56, 255), (0, 0, 85, 255), (0, 0, 113,255), (0, 0, 142, 255), (0, 0, 170, 255), (0, 0, 199, 255), (0, 0, 227, 255), (0, 0, 255, 255), (0, 28, 0, 255), (0, 28, 28, 255), (0, 28, 56, 255), (0, 28, 85, 255), (0, 28, 113,
255), (0, 28, 142, 255), (0, 28, 170, 255), (0, 28, 199, 255), (0, 28, 227, 255), (0, 28, 255, 255), (0, 56, 0, 255), (0, 56, 28, 255), (0,56, 56, 255), (0, 56, 85, 255), (0, 56, 113, 255), (0, 56, 142, 255), (0, 56, 170, 255), (0, 56, 199, 255), (0, 56, 227,
255), (0, 56, 255, 255), (0, 85, 0, 255), (0, 85, 28, 255), (0, 85, 56, 255), (0, 85, 85, 255), (0, 85, 113, 255), (0,85, 142, 255), (0, 85, 170, 255), (0, 85, 199, 255), (0, 85, 227, 255),(0,85,255, 255), (0, 113, 0, 255), (0, 113, 28, 255), (0, 113, 56,
255), (0, 113,85,255), (0, 113, 113, 255), (0, 113, 142, 255), (0, 113, 170, 255), (0, 113, 199,255), (0, 113, 227, 255), (0, 113, 255, 255), (0, 142, 0, 255), (0, 142, 28, 255), (0, 142, 56, 255), (0, 142, 85, 255), (0, 142, 113, 255), (0, 142, 142, 255),
(0, 142, 170, 255), (0, 142, 199, 255), (0, 142, 227, 255), (0, 142, 255, 255), (0, 170, 0, 255), (0, 170, 28, 255), (0, 170, 56, 255), (0, 170, 85, 255), (0,170, 113, 255), (0, 170, 142, 255), (0, 170, 170, 255), (0, 170, 199, 255), (0,170, 227, 255), (0,
170, 255, 255), (0, 199, 0, 255), (0, 199, 28, 255), (0, 199, 56, 255), (0, 199, 85, 255), (0, 199, 113, 255), (0, 199, 142, 255), (0, 199,170, 255), (0, 199, 199, 255), (0, 199, 227, 255), (0, 199, 255, 255), (0, 227,0, 255), (0, 227, 28, 255), (0, 227, 56,
255), (0, 227, 85, 255), (0, 227, 113,255), (0, 227, 142, 255), (0, 227, 170, 255), (0, 227, 199, 255), (0, 227, 227,255), (0, 227, 255, 255), (0, 255, 0, 255), (0, 255, 28, 255), (0, 255, 56, 255), (0, 255, 85, 255), (0, 255, 113, 255), (0, 255, 142, 255),
(0, 255, 170, 255), (0, 255, 199, 255), (0, 255, 227, 255), (0, 255, 255, 255), (28, 0, 0, 255), (28, 0, 28, 255), (28, 0, 56, 255), (28, 0, 85, 255), (28, 0, 113, 255), (28, 0, 142, 255), (28, 0, 170, 255), (28, 0, 199, 255), (28, 0, 227, 255), (28, 0, 255,255),
(28, 28, 0, 255), (28, 28, 28, 255), (28, 28, 56, 255), (28, 28, 85, 255), (28, 28, 113, 255), (28, 28, 142, 255), (28, 28, 170, 255), (28, 28, 199, 255), (28, 28, 227, 255), (28, 28, 255, 255), (28, 56, 0, 255), (28, 56, 28, 255), (28, 56, 56, 255), (28,
56, 85, 255), (28, 56, 113, 255), (28, 56, 142, 255), (28, 56, 170, 255), (28, 56, 199, 255), (28, 56, 227, 255), (28, 56, 255, 255), (28, 85, 0, 255), (28, 85, 28, 255), (28, 85, 56, 255), (28, 85, 85, 255), (28, 85,113, 255), (28, 85, 142, 255), (28, 85,
170, 255), (28, 85, 199, 255), (28, 85,227, 255), (28, 85, 255, 255), (28, 113, 0, 255), (28, 113, 28, 255), (28, 113,56, 255), (28, 113, 85, 255), (28, 113, 113, 255), (28, 113, 142, 255), (28, 113, 170, 255), (28, 113, 199, 255), (28, 113, 227, 255), (28,
113, 255, 255), (28, 142, 0, 255), (28, 142, 28, 255), (28, 142, 56, 255), (28, 142, 85, 255), (28,142, 113, 255), (28, 142, 142, 255), (28, 142, 170, 255), (28, 142, 199, 255), (28, 142, 227, 255), (28, 142, 255, 255), (28, 170, 0, 255), (28, 170, 28, 255),(28,
170, 56, 255), (28, 170, 85, 255), (28, 170, 113, 255), (28, 170, 142, 255), (28, 170, 170, 255), (28, 170, 199, 255), (28, 170, 227, 255), (28, 170, 255,255), (28, 199, 0, 255), (28, 199, 28, 255), (28, 199, 56, 255), (28, 199, 85, 255), (28, 199, 113, 255),
(28, 199, 142, 255), (28, 199, 170, 255), (28, 199, 199, 255), (28, 199, 227, 255), (28, 199, 255, 255), (28, 227, 0, 255), (28, 227,28, 255), (28, 227, 56, 255), (28, 227, 85, 255), (28, 227, 113, 255), (28, 227, 142, 255), (28, 227, 170, 255), (28, 227,
199, 255), (28, 227, 227, 255), (28,227, 255, 255), (28, 255, 0, 255), (28, 255, 28, 255), (28, 255, 56, 255), (28,255, 85, 255), (28, 255, 113, 255), (28, 255, 142, 255), (28, 255, 170, 255), (28, 255, 199, 255), (28, 255, 227, 255), (28, 255, 255, 255),
(56, 0, 0, 255), (56, 0, 28, 255), (56, 0, 56, 255), (56, 0, 85, 255), (56, 0, 113, 255), (56, 0, 142, 255), (56, 0, 170, 255), (56, 0, 199, 255), (56, 0, 227, 255), (56, 0, 255,255), (56, 28, 0, 255), (56, 28, 28, 255), (56, 28, 56, 255), (56, 28, 85, 255),
(56, 28, 113, 255), (56, 28, 142, 255), (56, 28, 170, 255), (56, 28, 199, 255), (56, 28, 227, 255), (56, 28, 255, 255), (56, 56, 0, 255), (56, 56, 28, 255), (56, 56, 56, 255), (56, 56, 85, 255), (56, 56, 113, 255), (56, 56, 142, 255), (56, 56, 170, 255), (56,
56, 199, 255), (56, 56, 227, 255), (56, 56, 255, 255), (56, 85, 0, 255), (56, 85, 28, 255), (56, 85, 56, 255), (56, 85, 85, 255), (56, 85,113, 255), (56, 85, 142, 255), (56, 85, 170, 255), (56, 85, 199, 255), (56, 85,227, 255), (56, 85, 255, 255), (56, 113,
0, 255), (56, 113, 28, 255), (56, 113,56, 255), (56, 113, 85, 255), (56, 113, 113, 255), (56, 113, 142, 255), (56, 113, 170, 255), (56, 113, 199, 255), (56, 113, 227, 255), (56, 113, 255, 255), (56, 142, 0, 255), (56, 142, 28, 255), (56, 142, 56, 255), (56,
142, 85, 255), (56,142, 113, 255), (56, 142, 142, 255)]

>>>

可以看到这里的数据是经过(原始数据/65535.0*255)调整过的。

这里在使用的时候可能比较方便,但是如果你是有读取过geotiff原始数据背景的,则需要小心。可能会有习惯性思维的陷阱。因为两者的类型都是short,但是实际的数值不一样(我就曾经以为gdal读取的是geotiff内部的原始ColorMap数据)。

3.略缩图集
http://www.gdal.org/gdal_datamodel.html
一个波段可能有0个或者更多的略缩图。每个略缩图被表示为一个“单独存在”的GDALRasterBand。略缩图大小将和处于下层的栅格不一样,但是在地理区域范围来说略缩图同样覆盖了整个区域。

略缩图主要用来在显示缩小分辨率时比读完整分辨率的图层数据更快。

波段会有一个HasArbitraryOverviews属性(如果图像在没有明显的略缩图图层的情况下,最高分辨率和低分辨率下都可以很有效率得进行读取,这个属性为true)这应用于某些fft编码图形或帮助那些在细微处理上需要效率的图形度过速度瓶颈的难关。(这句很别扭)

有两个地方可以建立金字塔
http://www.gdal.org/classGDALRasterBand.html#4dafcc399c9b4023f34cfdece3c6c574
函数原型

CPLErr GDALRasterBand::BuildOverviews      (       const char *        pszResampling,

        int      nOverviews,

        int *      panOverviewList,

        GDALProgressFunc      pfnProgress,

        void *      pProgressData

    )      [virtual]

用于建立波段略缩图。

如果制定数据集不支持这个操作,函数会返回CE_Failure。CPLGetLastErrorNo()将会返回CPLE_NotSupported值。

警告:对于单一波段的TIFF格式,不可能建立略缩图集,而且这样这个方法对于TIFF格式以及任何用在TIFF格式中建立的默认略缩图集的格式来说不都会正常工作。但是你可以在整个数据集中建立略缩图集(使用GDALDataset::BuildOverviews())。从实践角度来说这个方法相当没有用。

参数:

    pszResampling     "NEAREST", 平均值 或者"MODE"中的一个,用来控制重采样方法。

    nOverviews     一共需要建立多少个略缩图

    panOverviewList     一个数组,用来表示建立略缩图时缩小的参数列表

    pfnProgress     一个用来记录进程执行进度的回调函数(大概可以用来画进度条),如果不要处理进度,设置null

    pProgressData     传递到进度处理函数中的应用数据.

返回:

    成功返回CE_None,如果这个函数不能正常工作,返回CE_Failure .

既然这个函数没有意义,我们就看GDALDataset的同名函数
http://www.gdal.org/classGDALDataset.html#4925c147a1a168d01b7221e24869d5c3
函数原型:

CPLErr GDALDataset::BuildOverviews      (       const char *        pszResampling,

        int      nOverviews,

        int *      panOverviewList,

        int      nListBands,

        int *      panBandList,

        GDALProgressFunc      pfnProgress,

        void *      pProgressData

    )  

用于建立波段略缩图。

如果制定数据集不支持这个操作,方法会返回CE_Failure。CPLGetLastErrorNo()将会返回CPLE_NotSupported值。

这个方法的作用和C函数GDALBuildOverviews()作用是一样的。

参数:

        pszResampling     "NEAREST", 平均值 或者"MODE"中的一个,用来控制重采样方法

        nOverviews     一共需要建立多少个略缩图

        panOverviewList     一个数组,用来表示建立略缩图时缩小的参数列表

        nBand     在下面的波段号列表中需要建立略缩图的波段数,如果是0则建立所有波段

        panBandList     要建立略缩图的波段号列表

        pfnProgress     一个用来记录进程执行进度的回调函数(大概可以用来画进度条),如果不要处理进度,设置null

        pProgressData     传递到进度处理函数中的应用数据

返回:

    成功返回CE_None,如果这个函数不能正常工作,返回CE_Failure .

举一个例子,在所有波段上建立缩小比率为2,4,8的略缩图,需要这么写

   int       anOverviewList[3] = { 2, 4, 8 };

   poDataset->BuildOverviews( "NEAREST", 3, anOverviewList, 0, NULL, 

                              GDALDummyProgress, NULL );

 

-------------------------------------

对于这个略缩图集来说,其实就是传说中的栅格金字塔。这个可以从qgis的源码中看出来。建立金字塔所用的就是这个函数,而且到qgis中,调用的函数名字就叫建立金字塔。    

其实看起来这个函数其实并不是有太多的用处。虽然金字塔是一个很诱人的字样,但是在我研究后发现这个函数建立金字塔和jai的scale方法有本质区别。 jai的scale只是设置一个标志,不是真的去建,用时极少。然而这个BuildOverviews是真建。可以用python试一试。

>>> import gdal

>>> ds = gdal.Open("f:/gisdata/gtif/sd3.tif")

>>> help(ds.BuildOverviews)

Help on method BuildOverviews in module gdal:

BuildOverviews(self, resampling='NEAREST', overviewlist=None, callback=None, callback_data=None) method of gdal.Dataset instance

>>> ovlist = [2,4,8]

>>> ds.BuildOverviews(overviewlist=ovlist)

0

>>> 

可以看见在同路径下多了一个sd3.tif.ovr文件,大小1010k,那个就是金字塔所在。

多建,超过金字塔封顶也可以:

>>> ovlist = [2,4,8,16,32,64,128,256,1024]

>>> ds.BuildOverviews(overviewlist=ovlist)

0

>>>

这时文件大小又有所增加,变成1411k

多一层金字塔就多一点大小,建过的不会重复建。

>>> ovlist = [2,4,8]

>>> ds.BuildOverviews(overviewlist=ovlist)

0

>>> 

再次运行,ovr文件还是那么大。没有变化

单波段好像的确不能建立金字塔。

>>> import gdal

>>> ds = gdal.Open("f:/gisdata/gtif/asp.tif")

>>> dir(ds)

['AddBand', 'AdviseRead', 'BuildOverviews', 'FlushCache', 'GetDescription', 'GetDriver', 'GetGCPCount', 'GetGCPProjection', 'GetGCPs', 'GetGeoTransform', 'GetMetadata', 'GetProjection', 'GetProjectionRef', 'GetRasterBand', 'GetSubDatasets', 'RasterCount', 'RasterXSize',
'RasterYSize', 'ReadAsArray', 'ReadRaster', 'RefreshBandInfo', 'SetDescription', 'SetGCPs', 'SetGeoTransform', 'SetMetadata', 'SetProjection', 'WriteRaster', '__del__', '__doc__', '__init__', '__module__', '_band', '_o']

>>> help(ds.BuildOverviews)

Help on method BuildOverviews in module gdal:

BuildOverviews(self, resampling='NEAREST', overviewlist=None, callback=None, callback_data=None) method of gdal.Dataset instance

>>> ovlist = [2,4,8]

>>> ds.BuildOverviews(overviewlist=ovlist)

3

>>> 

返回是3,CE_Failure,失败。而且我方法是从dataset调的,不是band,这样也不行!!!

看来这个方法真的是没有什么用。除非真的需要在硬盘上实际建立金字塔(带实际数据的),且真的不心疼硬盘以及建立的时间。要不然还是用rasterio直接读比较好。
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: