# pythongdal教程之：地图代数与栅格数据的写入

ndvi=(nir-red)/(nir+red)

1. 将波段3读入数组data3，将波段2读入数组data2

2. 计算公式为：

mask = numeric.greater(data3 + data2, 0)

ndvi = numeric.choose(mask, (-99, (data3 – data2) / (data3 + data2)))

driver = indataset.getdriver()

create(, , , [], [])

outdataset = driver.create(filename, cols, rows, 1, gdt_float32)

outband = outdataset.getrasterband(1)

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)

nd = outband.getnodatavalue()

outband.flushcache()

outband.getstatistics(0, 1)

geotransform = indataset.getgeotransform()

outdataset.setgeotransform(geotransform )

proj = indataset.getprojection()

outdataset.setprojection(proj)

gdal.setconfigoption(‘hfa_use_rrd’, ‘yes’)

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)

6. 对输出图像：计算统计量，设定geotransform ：[minx, pixelwidth, 0, maxy, 0, pixelheight]，设定projection，建立pyramids

Posted in 未分类