Python 介绍三 ——Gdal

 

一、安装

安装GDAL可不是那么幸运的事儿，还是有点麻烦。Python2.3，python2.4的版本有安装包，而python2.5版本就没有安装包，得自己编译，因为gdal是C++写的模块，编译时必需安装vs.net2005，还是有点麻烦的。

一般来讲，必须包括三个文件/文件夹，gdal-1.4.0（是编译好的gdal）、gdal（这是一个路径的链接，将其链接到gdal-1.4.0文件夹中，只是一个简单的文本）拷到C:/Python25/Lib/site-packages目录下，注意gdal-1.4.0名字不能改（打开gdal看看，你就知道为什么），将gdal-140（gdal-140可从网上直接下的，可用于C++接口的，用于python时还必须与gdal-1.4.0的配合才行。）拷到一定的目录下，并将这一目录添加到PATH中既可。此外，如果你用的是python2.5，还必须安装numeric用于2.5才行，本来Numeric是用于2.4的，2.5得自已搞一个才行，否则装不上去。

二、概述

包括gdal和ogr的库，gdal支持超多的栅格的格式，因此很受欢迎的，如在设计数字图像处理的计算算法时很管用，ogr是支持矢量格式，ogr是建立在gdal的基础上，只有装了gdal才能用的。

三、功能介绍

先来个例子说明吧。

import gdal

dataset=gdal.Open( "a.jpg" )

band = dataset.GetRasterBand( 1 )

a=band.ReadAsArray( 100 , 100 , 5 , 5 , 10 , 10 )

Open>>gdal.Open()打开一个文件。

GetRasterBand(1)>>获得第一波段的数据。可见GDAL与地理数据关系很紧密的，也是按波段的数据来组织的。

ReadAsArray( 100 , 100 , 5 , 5 , 10 , 10 )>> 将读出图片某一波段一定范围的数字表示。

 

    再看看另一个程序，不完整，但足已说明其主要功能了。

Band1=self.dataset1.GetRasterBand(BandNumber1)

Band2=self.dataset2.GetRasterBand(BandNumber2)

Array1=Band1.ReadAsArray(0,0,self.Width,self.Height)

Array2=Band2.ReadAsArray(0,0,self.Width,self.Height)

Array3=Array1+Array2

driver=gdal.GetDriverByName("GTiff")             NewImage=driver.Create(self.OutFile,self.Wid,self.Heig,1,options=["INTERLEAVE=PIXEL"])            NewImage.WriteRaster(0,0,self.Wid,self.Heig,Array3.tostring(),self.Wid,self.Heig,band_list=[1])

 

    主要完成两个波段的相加。

GetDriverByName>>产生一个驱动，这是生成数据集的前提条件。

Create>>产生一个空的数据集。

WriteRaster>>将数据（Array3）写进数据集。

 

    再写个例子，涉及栅格数据的投影问题。

dataset=gdal.Open( "utm.tif" )

x=dataset.GetProjectionRef()

NewImage=driver.Create( "outpututm.tif" ,Width,Height, 3 ,options=[ "INTERLEAVE=PIXEL" ])

NewImage.WriteRaster( 0 , 0 ,Width,Height,Array1.tostring(),Width,Height,band_list=[ 1 ])

NewImage.SetGeoTransform( [ 592485 , 200 , 0 , 3494382 , 0 , - 200 ] )

NewImage.SetProjection(x)

以上过程，生成一幅图-- outpututm.tif，并且将它的投影设为与utm.tif的投影一样，其中， [ 592485 , 200 , 0 , 3494382 , 0 , - 200 ]是配准的结果，分别包括左上角的坐标，及X轴、y轴的分辨率。（这个与jpw,tfw,bpw的文件内容是一致的。）

GetProjectionRef()>>获得utm.tif图像的投影信息。

SetGeoTransform()>>设置左上角起始点的坐标，实际上就是配准的过程了。

SetProjection(x)>>设置投影。

不过，还是有点晕，这个投影的转换过程实际是没有涉及重采样的过程，也不知误差会有多大，还没验证。