以计算ndvi为例：
ndvi=(nir-red)/(nir+red)
其中nir为波段3，red为波段2
编程要点如下：
1. 将波段3读入数组data3，将波段2读入数组data2
2. 计算公式为：
3. 当data3和data2均为0时（例如用0表示nodata），会出现被0除的错误，导致程序崩溃。需要用mask配合choose将0值去掉
代码如下，仅有4行
data2 = band2.readasarray(0, 0, cols,rows).astype(numeric.float16)
data3 = band3.readasarray(0, 0, cols,rows).astype(numeric.float16)
mask = numeric.greater(data3 + data2, 0)
ndvi = numeric.choose(mask, (-99, (data3 – data2) / (data3 + data2)))
新建栅格数据集
将刚才计算得到的数据写入新的栅格数据集之中
首先要复制一份数据驱动：
driver = indataset.getdriver()
之后新建数据集
create(, , , [], [])
其中bands的默认值为1，gdaldatatype的默认类型为gdt_byte，例如
outdataset = driver.create(filename, cols, rows, 1, gdt_float32)
在这条语句的执行过程中，存储空间已经被分配到硬盘上了
在写入之前，还需要先引入波段对象
outband = outdataset.getrasterband(1)
波段对象支持直接写入矩阵，两个参数分别为x向偏移和y向偏移
outband.writearray(ndvi, 0, 0)
下面的例子总结了本次和上次的逐块写入方法
xblocksize = 64
yblocksize = 64
for i in range(0, rows, yblocksize):
if i + yblocksize < rows:
numrows = yblocksize
else:
numrows = rowsnumrows = rows –– ii
for j in range(0, cols, xblocksize):
if j + xblocksize < cols:
numcols = xblocksize
else:
numcols = cols – j
data = band.readasarray(j, i, numcols, numrows)
# do calculations here to create outdata array
outband.writearray(outdata, j, i)
band对象可以设定nodata值
outband.setnodatavalue(-99)
还可以读取nodata值
nd = outband.getnodatavalue()
计算band的统计量
首先用flushcache()把缓存数据写入磁盘
之后用getstatistics(, )计算统计量。如果approx_ok=1那么计算是基于pyramid的，如果force=0那么当整幅图都要被重读一遍的时候就不计算统计量了。
outband.flushcache()
outband.getstatistics(0, 1)
设定新图的地理参考点
如果新图与另一张图的地理参考信息完全一致，那就很简单了
geotransform = indataset.getgeotransform()
outdataset.setgeotransform(geotransform )
proj = indataset.getprojection()
outdataset.setprojection(proj)
建立pyramids
设定imagine风格的pyramids
gdal.setconfigoption(‘hfa_use_rrd’, ‘yes’)
强制建立pyramids
outdataset.buildoverviews(overviewlist=[2,4, 8,16,32,64,128])
图像的拼接
1. 对每张图：读取行数和列数，原点(minx,maxy)，像素长，像素宽，并计算坐标范围
maxx1 = minx1 + (cols1 * pixelwidth)
miny1 = maxy1 + (rows1 * pixelheight)
2. 计算输出图像的坐标范围：
minx = min(minx1, minx2, …) maxx = max(maxx1, maxx2, …)
miny = min(miny1, miny2, …) maxy = max(maxy1, maxy2, …)
3. 计算输出图像的行数和列数：
cols = int((maxx – minx) / pixelwidth)
rows = int((maxy – miny) / abs(pixelheight)
4. 建立并初始化输出图像
5. 对每张待拼接的图：计算offset值
xoffset1 = int((minx1 – minx) / pixelwidth)
yoffset1 = int((maxy1 – maxy) / pixelheight)
读入数据并按照上面计算的offset写入
6. 对输出图像：计算统计量，设定geotransform ：[minx, pixelwidth, 0, maxy, 0, pixelheight]，设定projection，建立pyramids
以上就是python gdal教程之：地图代数与栅格数据的写入的内容，更多相关内容请关注php中文网（www.php1.cn）！