5. GDAL python教程（4）——用GDAL读取栅格数据

GDAL原生支持超过100种栅格数据类型，涵盖所有主流GIS与RS数据格式，包括

  • 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

5.1. 导入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('HFA')
driver.Register()

打开已有的栅格数据集：

fn = 'aster.img'
ds = gdal.Open(fn, GA_ReadOnly)
if ds is None:
   print 'Could not open ' + 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(<index>)，其参数为波段的索引号 然后用ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>)，读出从(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 < rows:
    numRows = yBSize
else:
    numRows = rows – i
    for j in range(0, cols, xBSize):
        if j + xBSize < cols:
            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