GDAL生成Erdas Imagine

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

导入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]
adfGeoTransform[1]
adfGeoTransform[2]
adfGeoTransform[3]
adfGeoTransform[4]
adfGeoTransform[5]
注意栅格数据集的坐标一般都是以左上角为基准的。
下面的例子是从一个栅格数据集中取出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 < 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