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] /* 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)的矩阵。如果将矩阵大小设为1×1，就是读取一个像素了。但是这一方法只能将读出的数据放到矩阵中，就算只读取一个像素也是一样。例如：
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格式则按64×64像素平铺。
如果一行是一个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
以上就是python gdal教程之：用gdal读取栅格数据的内容，更多相关内容请关注php中文网（www.php1.cn）！