GDAL生成Erdas Imagine

时间:2022-11-06 11:21:32

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