GDAL安装及介绍

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

先看看这段[http://www.gdal.org/ GDAL主页]上的英文介绍吧！

简单地说，GDAL是一个操作各种栅格地理数据格式的库。包括读取、写入、转换、处理各种栅格数据格式（有些特定的格式对一些操作如写入等不支持）。它使用了一个单一的抽象数据模型就支持了大多数的栅格数据（GIS对栅格，矢量，3D数据模型的抽象能力实在令人叹服）。当然除了栅格操作，这个库还同时包括了操作矢量数据的另一个有名的库[http://ogr.maptools.org/ ogr]（ogr这个库另外介绍），这样这个库就同时具备了操作栅格和矢量数据的能力，买一送一，这么合算的买卖为什么不做^^。最最最重要的是这个库是跨平台的，开源的！如今这个库对各种数据格式的支持强大到令人啧啧的地步了。如果你对他的强大有什么怀疑的话，看看[http://www.gdal.org/formats_list.html 这里]一大串的GDAL所支持格式清单，吓到了吧！再看看它的主页最后那些使用了它作为底层数据处理的软件列表吧！其中你可以不知道GRASS，你也可以不知道Quantum GIS (QGIS)，但是你总该知道 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 Files\Microsoft Visual Studio .NET 2003\Vc7\bin\vcvars32.bat" 注册vc的编译环境，或者vs2010tools里面的命令行工具

打 开gdal文件夹下的nmake.opt修改GDAL_HOME = "C:\warmerda\bld"把路径改到需要把gdal安装的地方。不改也可以。这里需要添加python支持，所以修改PY_INST_DIR = $(GDAL_HOME)\pymod把路径改成python下的Lib\site-packages文件夹下。PYDIR   =    "C:\Software\Python24" 改成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在.net2005下只能顺利编译1.2，1.3以上版本不能顺利编译，有一个地方指针转换出错。可能是2005的编译器比以往的严厉一点吧。

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

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

 

= 快速开始 =
其实在主站的[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']
>>>
}}}
这样我们就打开了这个文件。并且我们可以看到可以供我们调用的函数们。现在我们不做修改，不做添加，所以只要带有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, 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)
>>>
}}}
我们就把图像中位于230,270，宽度10高度10的数据读取出来了。

我们看完了数据集的主要函数。似乎已经够用了。的确，如果只是为了显示图像，这些的确已经够了。但是如果需要更多信息，我们就不得不进入波段操作数据(实际上我们大多数时候都需要进入band获取信息)。下面我们现在来看看刚才读取出来的那个band有些什么东西可以供我们操作的。

{{{
>>> 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算出的最大最小值，是排除了无意义值后计算出来的最大最小值。

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