GDAL库介绍

5.1. GDAL库介绍#

  • 可能你不玩GIS，不懂这个库到底有什么用，或者和python有什么关系。但是你要玩GIS，RS，你就应当知道这个库的价值。就算你不玩GIS，我想这个库对你也应该有致命的吸引力。为什么？看下面的介绍吧！

  • 先看看这段[[http://www.gdal.org/|GDAL主页]]上的英文介绍吧！ {{http://www.gdal.org/gdalicon.png}} is a translator library for raster geospatial data formats that is released under an[[http://www.gdal.org/faq.html#license|X/MIT]]style[[http://www.opensource.org/|Open Source]]license by the[[http://www.osgeo.org/|Open Source Geospatial Foundation]]. As a library, it presents a[[http://www.gdal.org/gdal_datamodel.html|single abstract data model]]to the calling application for all supported formats. It also comes with a variety of useful[[http://www.gdal.org/gdal_utilities.html|commandline utilities]]for data translation and processing.

  • 简单地说，GDAL是一个操作各种栅格地理数据格式的库。包括读取、写入、转换、处理各种栅格数据格式（有些特定的格式对一些操作如写入等不支持）。它使用了一个单一的抽象数据模型就支持了大多数的栅格数据（GIS对栅格，矢量，3D数据模型的抽象能力实在令人叹服）。当然除了栅格操作，这个库还同时包括了操作矢量数据的另一个有名的库[[http://ogr.maptools.org/|ogr]]（ogr这个库另外介绍），这样这个库就同时具备了操作栅格和矢量数据的能力，买一送一，这么合算的买卖为什么不做^^。

  • 最最最重要的是这个库是跨平台的，开源的！如今这个库对各种数据格式的支持强大到令人啧啧的地步了。如果你对他的强大有什么怀疑的话，看看[[http://www.gdal.org/formats_list.html|这里]]一大串的GDAL所支持格式清单，吓到了吧！再看看它的主页最后那些使用了它作为底层数据处理的软件列表吧！其中你可以不知道[[http://grass.itc.it/index.html|GRASS]]，你也可以不知道[[http://qgis.sourceforge.net/|Quantum GIS (QGIS)]]，但是你总该知道[[http://earth.google.com/|Google Earth]]吧！不知道？赶快下一个去玩玩－－会当临绝顶，一览众山小！

  • 有人说我又不玩GIS。不错，但是，你即使不玩GIS，这个库也是满有用的。首先，哪个库支持这么多栅格（图片）格式，哪个库在C/C++/python/ruby/VB/java/C#(这个暂时不完全支持)下都能用，而且都一样用？退一步讲，3S软件又不一定要用在3S下（很多医学影像就是用PCI软件来处理的）。再退一步，你的生活即使和3S一点关系都没有，栅格数据又不单单只有GIS下才用到。你大可用这个库来读取jpg，gif，tif，xpm等格式。而且对各种格式支持得不是一般的好，很大一部分非标准格式照样支持得非常好。我曾经在java下玩过jai，以及一系列jai的扩展库，一些图像格式在很多图片浏览器中都可以正确读取（有的甚至不是非标准格式），用jai死活就读不出来！

  • 这个库的python版和其他的python库结合的很好。最直接、明显的支持是使用Numeric库来进行数据读取和操作。各种矩阵魔术可以发挥得淋漓尽致（图像其实就是矩阵）。而且按我的观点，python对矩阵的操作比其他的语言有明显的优势。写出来的东西比其他语言写出来的短小的多，而且好看得多。并且python的弱类型在处理栅格数据格式类型的时候代码量比强类型的语言少了数倍（不用double，byte，short等等分开处理，这简直就是先天上的优势）。所以我就喜欢用python做图像的处理。所以就连GIS界的微软ESRI也直接在ARCGIS9中用python来作栅格数据的导入导出。一句话，真是太方便啦！

= 安装 = == windows下的安装 ==

  • 官方安装文档在[[http://www.remotesensing.org/gdal/gdal_building.html|这里]]。下面是我自己的实践步骤：

  • 先去http://www.gdal.org/dl/下一个版本，解压。 打开控制台，输入： “D:Program FilesMicrosoft Visual Studio .NET 2003Vc7binvcvars32.bat” 注册vc的编译环境。

  • 打 开gdal文件夹下的nmake.opt修改GDAL_HOME = “C:warmerdabld”把路径改到需要把gdal安装的地方。不改也可以。这里需要添加python支持，所以修改PY_INST_DIR = $(GDAL_HOME)pymod把路径改成python下的Libsite-packages文件夹下。PYDIR = “C:SoftwarePython24” 改成python的安装路径。 下面的参数爱改什么就把前面的#删除（要看您有没有那些库的源码），注意一下路径就可以了。我是都没改。 后面就依次运行

{{{ nmake /f makefile.vc nmake /f makefile.vc install nmake /f makefile.vc devinstall }}} 最后最后，还要去GDAL_HOME目录下的bin文件夹下把gdal13.dll （也有可能是gdal12.dll）copy到PY_INST_DIR路径下

  • 到此处就完成安装gdal（python）的工作。

  • 最后需要`注意一下`，gdal在vc++.net2005下只能顺利编译1.2，1.3系列的版本不能顺利编译，有一个地方指针转换出错。可能是2005的编译器比以往的严厉一点吧。但是vc++.net2005却可以正常编译1.4版本，不过要支持Python，需要改一个地方：

  • 打开Pymod目录下的makefile.vc，然后找到link这两行：

{{{
link /dll /def:_gdal.def $(OBJ) ../gdal_i.lib /LIBPATH(PYDIR)/libs

/out(PYGDAL_DLL)

5.1.1. }}}#

  • 然后在下面添加一行

{{{

if exist $(PYGDAL_DLL).manifest mt -manifest $(PYGDAL_DLL).manifest -outputresource:$(PYGDAL_DLL);2

5.1.2. }}}#

不加这行，会出现一个“找不到MSVCR80.DLL”的错误，这个错误是由于VS.net2005强制进行manifest验证造成的。而且这个manifest的问题相当棘手，对发布也有极坏的影响。

另外,安装了QGIS,对编译也有一些影响,主要是proj库的冲突,导致一个找不到”d:/program.obj”文件的错误,如果你有静态编译过proj,那么你可以打开nmake.opt修改有关proj的设置,如果搞不定,就卸载QGIS,然后编译,编译后再安装QGIS.呵呵,还好QGIS的体积没有ArcGIS那么可怕.

== linux下的安装 ==

  • linux下的安装就更简单了。直接

{{{ ./configure make su make install ldconfig }}} 就ok(默认就已经支持python)。当然在第一步的时候需要看看是否依赖的库都安装了。如果缺少，就去安装一个。如果对configure的条件不理解，就用./configure –help看看具体情况。

== 安装其他驱动 ==

  • 这里讲一个安装hdf4的驱动的例子（默认情况下gdal是不安装hdf4的），其他驱动应该和这个也差不了多少吧，可以作为其他的参考。完整步骤如下：

  • 在windows下的安装：

从[[ftp://ftp.ncsa.uiuc.edu/HDF/HDF/HDF_Current/bin/windows/|ftp://ftp.ncsa.uiuc.edu/HDF/HDF/HDF%5FCurrent/bin/windows/]]下载42r1-win.ZIP，解压。

  • 编辑gdal根目录下的nmake.opt，找到“# Uncomment the following and update to enable NCSA HDF Release 4 support.”这一行

  • 把下面两行前面的#去掉，然后改成：

{{{ HDF4_DIR = D:warmerda42r1-winrelease #HDF4_LIB = /LIBPATH:$(HDF4_DIR)lib hd421m.lib HDF4_LIB = $(HDF4_DIR)dllhd421m.lib $(HDF4_DIR)dllhm421m.lib

. $(HDF4_DIR)libhd421.lib $(HDF4_DIR)libhm421.lib

}}}

用HDF4_LIB=/LIBPATH：这种形式似乎可以建立gdal的库，但是往下编译会出错。而且要把$(HDF4_DIR)dll和$(HDF4_DIR)lib拷贝到同一个目录下，不然会提示找不到库

  • 你也可以试一试在`D:Program FilesMicrosoft Visual Studio .NET 2003Common7Toolsvsvars32.bat`文件中添加HDF4_LIB路径到“@set LIB=”这行的末尾（不要忘记;的分割符）。

  • 然后找一下”INC=”这行,把 -I$(HDF4_DIR)include加到下一行的末尾（应该也可以在vsvars32.bat中添加路径，不过要重启命令行）。

然后编译吧！祝你好运。

  • `注意：`上面的HDF4_DIR 是我本机的路径，你要根据你自己的路径进行设置（想起我的一个老师说过的话：“抄人家的作业可以，不要连名字也一起抄走啊 ” :) ），下面的$(HDF4_DIR)可以不用改，那个是变量，会自动替代HDF4_DIR 路径。

编译成功后，要HDF4能运行，还需要两个库，一个是zlib,一个是szip，可以到下面两个链接去下载一个

把这两个库下载后解压，然后设置PATH系统变量，使得它们在默认状态下也可以被动态链接成功 。

  • 在Linux下比在Windows下简单(好象永远如此;-)：

  • 只要用./configure –help察看一下打开HDF4的编译开关（包括库路径，头文件路径等，看清楚），然后在./configure 后面加上那个开关以及hdf4的安装路径后就可以了。在configure后gdal会提示是否支持HDF4。 编译后也要把zlib和szip的动态链接库设置好 。

  • 到此你已经可以用C/C++来操作gdal读写hdf4的格式了！

  • 最后，为了让Python能够支持hdf的读写，别忘了把重新生成安装后的pymod目录下的内容，还有gdal13，还有那两个hdf的库，还有zlib，szip的库拷贝到Python的Libsite-packages目录下 。

== 下载 ==

  • 如果你实在玩不转，可以在这里下载已经编译好的[[attachment:bld.zip|gdal1.3.2程序库]] 以及其[[attachment:extlib.zip|依赖的其他库]]，其中包括hdf4，hdf5支持，以及proj，geos插件。注意，这里的geos是静态链接的，注意版权（geos是LGPL的license）。hdf4和hdf5用的是release版本。[[attachment:nmake.opt|这里]]是我的nmake配置文件，你可以对照你的实际情况参考一下。

= 快速开始 =

  • 其实在主站的[[http://www.gdal.org/gdal_tutorial.html|教程]]里已经有python的示例了。但是我们还是按照自己的思路来开始吧。

  • 第一步就是打开一个数据集。对于“数据集”这个名词大家可能不会太习惯，但是对于一般的格式来说，一个“数据集”就是一个文件，比如一个gif文件就是一个以gif为扩展名的文件。但是对于众多RS数据来说，一个数据集包含的绝对不仅仅是一个文件。对于很多RS数据，他们把一张图像分成数个图像文件，然后放在一个文件夹中，用一些额外的文件来组织它们之间的关系，形成一个“数据集”。如果你不理解，那么就算了，当成jpg或者gif文件好了。

  • 下面我们打开一个tiff文件(GeoTIFF)。这个文件是我从GRASS的示例数据spearfish中导出的一个同名影像数据。

{{{ >>> import gdal >>> dataset = gdal.Open(“j:/gisdata/gtif/spot.tif”) >>> dir(dataset) [‘AddBand’, ‘AdviseRead’, ‘BuildOverviews’, ‘FlushCache’, ‘GetDescription’, ‘Get Driver’, ‘GetGCPCount’, ‘GetGCPProjection’, ‘GetGCPs’, ‘GetGeoTransform’, ‘GetMe tadata’, ‘GetProjection’, ‘GetProjectionRef’, ‘GetRasterBand’, ‘GetSubDatasets’,

‘RasterCount’, ‘RasterXSize’, ‘RasterYSize’, ‘ReadAsArray’, ‘ReadRaster’, ‘Refr

eshBandInfo’, ‘SetDescription’, ‘SetGCPs’, ‘SetGeoTransform’, ‘SetMetadata’, ‘Se tProjection’, ‘WriteRaster’, ‘__del__’, ‘__doc__’, ‘__init__’, ‘__module__’, ‘_b and’, ‘_o’] >>> }}}

  • 这样我们就打开了这个文件。并且我们可以看到可以供我们调用的函数们(更具体的API列表可以看[[http://www.gdal.org/classGDALDataset.html|这里]])。现在我们不做修改，不做添加，所以只要带有Set开头的函数以及有Write开头的函数我们暂时都不管。因为RS影像必然要和地理上的位置挂上钩，才能把图像正确铺展到一个坐标系中。其中的信息和对应关系有点复杂，不适合在快速开始中介绍，我们暂时也先不管。这里需要注意的就是几个函数。

  • !GetDescription 获得栅格的描述信息。

{{{ >>> dataset.GetDescription() ‘j:/gisdata/gtif/spot.tif’ >>> }}}

  • 看来这里的图像描述是图像的路径名，但是这是和各种不同数据集相关的，不同数据集可能有不同的描述。这要看读取驱动的实现作者的高兴了。

  • !RasterCount 获得栅格数据集的波段数。

  • !GetRasterBand 获得栅格数据集的波段。

{{{ >>> dataset.RasterCount 1 >>> band = dataset.GetRasterBand(1) >>> }}}

  • 这里需要解释的是Band这个词。这个词可以翻译成“波段”，“通道”等等。我这里把它统一称为“波段”。因为遥感卫星的传感器有很多个。一个传感器只负责接收一个频率范围的地物反射光波，一个频率范围的光波记录称为一个波段。是不是晕了？其实说得简单一点。其实你可以把波段看成红绿蓝几种颜色。图像不是分RGB三色吗？把R，G，B值都提取出来成为三个表。R值表就是波段一，G值表就是波段二，B值表就是波段三。

  • 这里我们看到这张图只有一个波段(一种颜色)。就可以把它看成是一个灰度图(类似黑白照片)。如果!RasterCount是3，就有可能是彩色图。如果!RasterCount是比3大的数，恭喜你，你看到一张遥感影像。有很多卫星的传感器大于3个，比如TM就有7个波段，不仅有可见光，还有红外等其他非可见光。，所以，波段一般比RGB能表达的丰富地多。不过这样一来就需要我们从中挑出3个波段然后组合成RGB，当然这样就有可能使图像显示出来的东西不像平常我们看到的那样。这样安排是因为对科学有帮助(一些波段在科学家眼里比真实的彩色照片有价值)。不理解就跳过，很正常，我第一次听这种东西也觉得很玄:)

  • 这里我们获取了第一个波段(红色值组成的表)。注意！这里的波段获取和通常的C数组获取不一样，开始是1不是0。获取了波段，我们就可以在下面的操作中读取这个波段的所有数值。

  • RasterXSize 图像的宽度(X方向上的像素个数)

  • RasterYSize 图像的高度(Y方向上的像素个数)

{{{ >>> dataset.RasterXSize 950 >>> dataset.RasterYSize 700 >>> }}} 可以看出我们的图像大小是950*700。还是很小的一张图。

  • !ReadRaster 读取图像数据(以二进制的形式)

  • !ReadAsArray 读取图像数据(以数组的形式)

{{{ >>> help(dataset.ReadRaster) Help on method ReadRaster in module gdal: ReadRaster(self, xoff, yoff, xsize, ysize, buf_xsize=None, buf_ysize=None, buf_t ype=None, band_list=None) method of gdal.Dataset instance >>> help(dataset.ReadAsArray) Help on method ReadAsArray in module gdal: ReadAsArray(self, xoff=0, yoff=0, xsize=None, ysize=None) method of gdal.Dataset

instance

>>>
}}}
 * 这两个函数很重要，它们直接读取图像的数据，可以看到两个函数的帮助中有一大溜的参数。解释一下： xoff,yoff,xsize,ysize 你可能不想读取整张图像。只想读取其中的一部分。那么就用xoff，yoff指定想要读取的部分原点位置在整张图像中距离全图原点的位置。用xsize和ysize指定要读取部分图像的矩形大小。
 * `buf_xsize buf_ysize` 你可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高，gdal将帮你缩放到这个大小。
 * `buf_type` 如果你要读取的图像的数据类型不是你想要的(比如原图数据类型是short，你要把它们缩小成byte)，就可以设置它。
 * `band_list *  这就适应上面多波段的情况。你可以指定读取的波段序列。要哪几个波段，不要哪几个波段，你说了算。

  • 举个例子吧：

{{{ >>> dataset.ReadAsArray(230,270,10,10) array([[255, 255, 255, 232, 232, 255, 255, 255, 255, 222],

[255, 255, 255, 255, 255, 255, 210, 110, 11, 122], [255, 255, 255, 255, 255, 255, 210, 255, 11, 243], [201, 255, 255, 255, 255, 200, 200, 110, 122, 243], [111, 211, 255, 201, 255, 255, 100, 11, 132, 243], [255, 100, 100, 100, 110, 100, 110, 111, 122, 243], [255, 255, 255, 255, 255, 255, 122, 222, 255, 255], [255, 255, 255, 255, 255, 255, 243, 243, 255, 255], [255, 255, 255, 255, 255, 255, 255, 255, 255, 255], [255, 255, 255, 255, 255, 255, 255, 255, 255, 255]],’b’)

>>> dataset.ReadRaster(230,270,10,10)
'\xff\xff\xff\xe8\xe8\xff\xff\xff\xff\xde\xff\xff\xff\xff\xff\xff\xd2n\x0bz\xff\
xff\xff\xff\xff\xff\xd2\xff\x0b\xf3\xc9\xff\xff\xff\xff\xc8\xc8nz\xf3o\xd3\xff\x
c9\xff\xffd\x0b\x84\xf3\xffdddndnoz\xf3\xff\xff\xff\xff\xff\xffz\xde\xff\xff\xff
\xff\xff\xff\xff\xff\xf3\xf3\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff
\xff\xff\xff\xff\xff\xff\xff\xff\xff'
>>>
}}}
 * 我们就把图像中位于230,270，宽度10高度10的数据读取出来了。
 * 我们看完了数据集的主要函数。似乎已经够用了。的确，如果只是为了显示图像，这些的确已经够了。但是如果需要更多信息，我们就不得不进入波段操作数据(实际上我们大多数时候都需要进入band获取信息)。下面我们现在来看看刚才读取出来的那个band有些什么东西可以供我们操作的(具体的API列表看[[http://www.gdal.org/classGDALRasterBand.html|这里]])。

{{{ >>> dir(band) [‘AdviseRead’, ‘Checksum’, ‘ComputeBandStats’, ‘ComputeRasterMinMax’, ‘DataType’ , ‘Fill’, ‘FlushCache’, ‘GetDefaultHistogram’, ‘GetDescription’, ‘GetHistogram’,

‘GetMaximum’, ‘GetMetadata’, ‘GetMinimum’, ‘GetNoDataValue’, ‘GetOffset’, ‘GetO

verview’, ‘GetOverviewCount’, ‘GetRasterColorInterpretation’, ‘GetRasterColorTab le’, ‘GetScale’, ‘GetStatistics’, ‘ReadAsArray’, ‘ReadRaster’, ‘SetDefaultHistog ram’, ‘SetDescription’, ‘SetMetadata’, ‘SetNoDataValue’, ‘SetRasterColorInterpre tation’, ‘SetRasterColorTable’, ‘WriteArray’, ‘WriteRaster’, ‘XSize’, ‘YSize’, ‘ __doc__’, ‘__init__’, ‘__module__’, ‘_o’] >>> }}}

  • 挑几个有用的吧。

{{{ >>> band.XSize 950 >>> band.YSize 700 >>> band.DataType 1 >>> }}}

  • 不用解释了吧，波段图像的宽和高(象元为单位)。!DataType，图像中实际数值的数据类型。具体数据类型定义在gdalconst模块里。使用的时候用import gdalconst引入。

{{{ >>> import gdalconst >>> dir(gdalconst) [‘CE_Debug’, ‘CE_Failure’, ‘CE_Fatal’, ‘CE_None’, ‘CE_Warning’, ‘CPLES_Backslash Quotable’, ‘CPLES_CSV’, ‘CPLES_SQL’, ‘CPLES_URL’, ‘CPLES_XML’, ‘CPLE_AppDefined’ , ‘CPLE_AssertionFailed’, ‘CPLE_FileIO’, ‘CPLE_IllegalArg’, ‘CPLE_NoWriteAccess’ , ‘CPLE_None’, ‘CPLE_NotSupported’, ‘CPLE_OpenFailed’, ‘CPLE_OutOfMemory’, ‘CPLE _UserInterrupt’, ‘CXT_Attribute’, ‘CXT_Comment’, ‘CXT_Element’, ‘CXT_Literal’, ‘ CXT_Text’, ‘DCAP_CREATE’, ‘DCAP_CREATECOPY’, ‘DMD_CREATIONDATATYPES’, ‘DMD_CREAT IONOPTIONLIST’, ‘DMD_EXTENSION’, ‘DMD_HELPTOPIC’, ‘DMD_LONGNAME’, ‘DMD_MIMETYPE’ , ‘GA_ReadOnly’, ‘GA_Update’, ‘GCI_AlphaBand’, ‘GCI_BlackBand’, ‘GCI_BlueBand’, ‘GCI_CyanBand’, ‘GCI_GrayIndex’, ‘GCI_GreenBand’, ‘GCI_HueBand’, ‘GCI_LightnessB and’, ‘GCI_MagentaBand’, ‘GCI_PaletteIndex’, ‘GCI_RedBand’, ‘GCI_SaturationBand’ , ‘GCI_Undefined’, ‘GCI_YellowBand’, ‘GDT_Byte’, ‘GDT_CFloat32’, ‘GDT_CFloat64’,

‘GDT_CInt16’, ‘GDT_CInt32’, ‘GDT_Float32’, ‘GDT_Float64’, ‘GDT_Int16’, ‘GDT_Int

32’, ‘GDT_TypeCount’, ‘GDT_UInt16’, ‘GDT_UInt32’, ‘GDT_Unknown’, ‘GF_Read’, ‘GF_ Write’, ‘GPI_CMYK’, ‘GPI_Gray’, ‘GPI_HLS’, ‘GPI_RGB’, ‘GRA_Bilinear’, ‘GRA_Cubic ‘, ‘GRA_CubicSpline’, ‘GRA_NearestNeighbour’, ‘__builtins__’, ‘__doc__’, ‘__file __’, ‘__name__’] >>> }}}

  • 那些GDT开头的就是数值数据类型。

{{{ >>> band.GetNoDataValue() 65535.0 >>> band.GetMaximum() >>> band.GetMinimum() >>> band.ComputeRasterMinMax() (1.0, 255.0) >>> }}}

  • Maximum 是表示在本波段数值中最大的值，Minimum当然就是表示本波段中最小的值啦。我们可以看到在一开始这两个都没有值。因为对于文件格式不会有固有的最大最小值。所以我们通过函数!ComputeRasterMinMax计算得到了。注意！这里的最大最小值不包括“无意义值”！也就是上面显示的!NoDataValue。需要解释一下“无意义值”。不要以为0或者255在任何情况下都无意义。在很多情况下0，255需要和其他值一样表示一个实际意义。虽然可能它最终会被显示得和黑色一样。而一些位置上的点要表示的意思是“什么也不是”，它在那个位置上只是为了占一个位置，使得整体图像看起来像个矩形而已。在做实际应用的时候两种值的处理将会完全不一样。所以需要设置无意义值，来和其他的值区别开来。而用!ComputeRasterMinMax算出的最大最小值，是排除了无意义值后计算出来的最大最小值。

{{{ >>> band.GetRasterColorInterpretation() 2 >>> gdalconst.GCI_PaletteIndex 2 >>> colormap = band.GetRasterColorTable() >>> dir(colormap) [‘Clone’, ‘GetColorEntry’, ‘GetColorEntryAsRGB’, ‘GetCount’, ‘GetPaletteInterpre tation’, ‘SetColorEntry’, ‘__del__’, ‘__doc__’, ‘__init__’, ‘__module__’, ‘__str __’, ‘_o’, ‘own_o’, ‘serialize’] >>> colormap.GetCount() 256 >>> colormap.GetPaletteInterpretation() 1 >>> gdalconst.GPI_RGB 1 >>> for i in range(colormap.GetCount()): … print colormap.GetColorEntry(i), … (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, 1 42, 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, 19 9, 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,

  1. (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, 2

55) (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,

  1. (28, 56, 255, 255) (28, 85, 0, 255) (28, 85, 28, 255) (28, 85, 56, 255) (2

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

  1. (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, 25 5) (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, 5 6, 255) (56, 28, 85, 255) (56, 28, 113, 255) (56, 28, 142, 255) (56, 28, 170, 25 5) (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, 1 42, 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, 25 5) (56, 142, 28, 255) (56, 142, 56, 255) (56, 142, 85, 255) (56, 142, 113, 255) (56, 142, 142, 255) >>> }}}

  • 通过`!GetRasterColorInterpretation`，我们知道我们的图像是一个颜色表索引的图像而不是纯粹的黑白灰度图像(!PaletteIndex，其他的颜色模型，可以察看gdalconst模块中GCI打头的枚举值)。这意味着我们读出的数据有可能不是真实的数据。这些数据只是一个个实际数据的索引。真实数据存储在另一个表中。我们通过!ReadRaster读出的数据值只是对应到这个表的一个索引而已。我们需要通过读出这些数据，并在真实数据表中找出真实数据，重新组织成一个RGB表才能用来绘制。如果我们不经过对应，我们绘制出来的东西可能什么东西都不是。

  • 用`!GetRasterColorTable`获得了颜色表，通过`!GetPaletteInterpretation`我们知道我们获得的颜色表是一个RGB颜色表。GDAL支持多种颜色表，具体可以参考gdalconst模块中GPI打头的枚举值。然后我们可以通过!GetCount获得颜色的数量。通过!GetColorEntry获得颜色表中的值。这里的颜色值都是一个4值的元组。里面有意义的只有前三个(如果颜色模型是GPI_RGB, GPI_HLS，都使用前3个，如果采用GPI_CMYK，则4个值都有意义了)。

{{{ >>> 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 >>> help(band.ReadRaster) Help on method ReadRaster in module gdal: ReadRaster(self, xoff, yoff, xsize, ysize, buf_xsize=None, buf_ysize=None, buf_t ype=None) method of gdal.Band instance >>> }}}

  • 显然，band里的`!ReadAsArray`参数显然比dataset里面的要好用，而`!ReadRaster`则差不多。但是!ReadAsArray读出的是数组，可以用Numeric模块进行矩阵魔法。!ReadRaster读出的是二进制，虽然可以直接绘制，但是对于一些绘图API来说，对[[RRR…][GGG…][BBB…]]表的处理明显不如[[RGB][RGB]…]，有的甚至不支持。虽然可以用struct.unpack来拆封，可效率上就差很多(而且拆封出来还是数组)。数组在矩阵魔法的控制之下则会显得十分方便快捷，最后用tostring直接转化称为二进制绘制，速度也相当快。

  • 好了，快速开始已经使我们可以初步看清楚了gdal中图像的组织。下面用一句话总结一下：波段组成图像，波段指挥颜色。

= 反馈 = 如果您发现我写的东西中有问题，或者您对我写的东西有意见，请登陆[[http://groups.google.com/group/geosings?hl=zh-CN|这个论坛]]

Contents