pythongdal教程之：用gdal读取栅格数据

ArcInfo grids, ArcSDE raster, Imagine, Idrisi, ENVI, GRASS, GeoTIFF

HDF4, HDF5

USGS DOQ, USGS DEM

ECW, MrSID

TIFF, JPEG, JPEG2000, PNG, GIF, BMP

完整的支持列表可以参考http://www.gdal.org/formats_list.html

导入GDAL支持库

旧版本（1.5以前）：import gdal, gdalconst

新版本（1.6以后）：from osgeo import gdal, gdalconst

gdal和gdalconst最好都要导入，其中gdalconst中的常量都加了前缀，力图与其他的module冲突最小。所以对gdalconst你可以直接这样导入：from osgeo.gdalconst import *

GDAL数据驱动，与OGR数据驱动类似，需要先创建某一类型的数据驱动，再创建响应的栅格数据集。

一次性注册所有的数据驱动，但是只能读不能写：gdal.AllRegister()

单独注册某一类型的数据驱动，这样的话可以读也可以写，可以新建数据集：

driver = gdal.GetDriverByName(&＃39;HFA&＃39;)

driver.Register()

打开已有的栅格数据集：

fn = &＃39;aster.img&＃39;

ds = gdal.Open(fn, GA_ReadOnly)

if ds is None:

print &＃39;Could not open &＃39; + fn

sys.exit(1)

读取栅格数据集的x方向像素数，y方向像素数，和波段数

cols = ds.RasterXSize

rows = ds.RasterYSize

bands = ds.RasterCount

注意后面没有括号，因为他们是属性(properties)不是方法(methods)

读取地理坐标参考信息(georeference info)

GeoTransform是一个list，存储着栅格数据集的地理坐标信息

adfGeoTransform[0] /* top left x 左上角x坐标*/

adfGeoTransform[1] /* w--e pixel resolution 东西方向上的像素分辨率*/

adfGeoTransform[2] /* rotation, 0 if image is "north up" 如果北边朝上，地图的旋转角度*/

adfGeoTransform[3] /* top left y 左上角y坐标*/

adfGeoTransform[4] /* rotation, 0 if image is "north up" 如果北边朝上，地图的旋转角度*/

adfGeoTransform[5] /* n-s pixel resolution 南北方向上的像素分辨率*/

注意栅格数据集的坐标一般都是以左上角为基准的。

下面的例子是从一个栅格数据集中取出Geotransform作为一个list，然后读取其中的数据

geotransform = ds.GetGeoTransform()

originX = geotransform[0]

originY = geotransform[3]originY = geotransform[3]

pixelWidth = geotransform[1]

pixelHeight = geotransform[5]

计算某一坐标对应像素的相对位置(pixel offset)，也就是该坐标与左上角的像素的相对位置，按像素数计算，计算公式如下：

xOffset = int((x – originX) / pixelWidth)

yOffset = int((y – originY) / pixelHeight)

读取某一像素点的值，需要分两步

首先读取一个波段(band)：GetRasterBand()，其参数为波段的索引号

然后用ReadAsArray(, , , )，读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。如果将矩阵大小设为1X1，就是读取一个像素了。但是这一方法只能将读出的数据放到矩阵中，就算只读取一个像素也是一样。例如：

band = ds.GetRasterBand(1)

data = band.ReadAsArray(xOffset, yOffset, 1, 1)

如果想一次读取一整张图，那么将offset都设定为0，size则设定为整个图幅的size，例如：

data = band.ReadAsArray(0, 0, cols, rows)

但是要注意，从data中读取某一像素的值，必须要用data[yoff, xoff]。注意不要搞反了。数学中的矩阵是[row,col]，而这里恰恰相反！这里面row对应y轴，col对应x轴。

注意在适当的时候释放内存，例如band = None 或者dataset = None。尤其当图很大的时候

如何更有效率的读取栅格数据？显然一个一个的读取效率非常低，将整个栅格数据集都塞进二维数组也不是个好办法，因为这样占的内存还是很多。更好的方法是按块(block)来存取数据，只把要用的那一块放进内存。本周的样例代码中有一个utils模块，可以读取block大小。

例如：

import utils

blockSize = utils.GetBlockSize(band)

xBlockSize = blockSize[0]

yBlockSize = blockSize[1]

平铺(tiled)，即栅格数据按block存储。有的格式，例如GeoTiff没有平铺，一行是一个block。Erdas imagine格式则按64x64像素平铺。

如果一行是一个block，那么按行读取是比较节省资源的。

如果是平铺的数据结构，那么设定ReadAsArray()的参数值，让它一次只读入一个block，就是效率最高的方法了。例如：

rows = 13, cols = 11, xBSize = 5, yBSize = 5

for i in range(0, rows, yBSize):

if i + yBSize

numRows = yBSize

else:

numRows = rows – i

for j in range(0, cols, xBSize):

if j + xBSize

numCols = xBSize

else:

numCols = colsnumCols = cols – j

data = band.ReadAsArray(j, i, numCols, numRows)

这一段代码具有通用性，可以时常拿来用的。

下面介绍一点二维数组的处理技巧

这里要用到两个库，Numeric和numpy。Numeric比较老了，FWTools用它。自己安装配置的话还是配功能更强的numpy。

数据类型转换：

data = band.ReadAsArray(j, i, nCols, nRows)

data = data.astype(Numeric.Float) # Numeric

data = data.astype(numpy.float) # numpy

或者简单点只写一句

data = band.ReadAsArray(j, i, nCols, nRows).astype(Numeric.Float)

掩膜mask

这是Numeric和numpy库的功能，输入一个数组和条件，输出一个二值数组。例如

mask = Numeric.greater(data, 0)mask = Numeric.greater(data, 0)

>>> a = Numeric.array([0, 4, 6, 0, 2])

>>> print a

[0 4 6 0 2]

>>> mask = Numeric.greater(a, 0)

>>> print mask

[0 1 1 0 1]

数组求和

>>> a = Numeric.array([0, 4, 6, 0, 2])

>>> print a>>> print a

[0 4 6 0 2]

>>> print Numeric.sum(a)

12

如果是二维数组，那sum就会返回一个一维数组

>>> b = Numeric.array([a, [5, 10, 0, 3, 0]])

>>> print b

[[ 0 4 6 0 2]

[ 5 10 0 3 0]]

>>> print Numeric.sum(b)>>> print Numeric.sum(b)

[ 5 14 6 3 2]

所以，二维数组的求和就要这样

>>> print Numeric.sum(Numeric.sum(b))

30

这里有一个小技巧，统计大于0的像素个数，可以联合运用mask和sum两个函数

>>> print a

[0 4 6 0 2]

>>> mask = Numeric.greater(a, 0)

>>> print mask

[0 1 1 0 1]

>>> print Numeric.sum(mask)

3

以上就是python gdal教程之：用gdal读取栅格数据的内容，更多相关内容请关注PHP中文网（www.php1.cn）！