# 6. GDAL python教程（5）——地图代数与栅格数据的写入¶

## 6.1. 以计算NDVI为例：¶

```NDVI=(NIR-RED)/(NIR+RED)
```

1. 将波段3读入数组data3，将波段2读入数组data2
2. 计算公式为： ndvi = (data3 – data2) / (data3 + data2)

```data2 = band2.ReadAsArray(0, 0, cols,rows).astype(Numeric.Float16)
mask = Numeric.greater(data3 + data2, 0)
ndvi = Numeric.choose(mask, (-99, (data3 - data2) / (data3 + data2)))
```

## 6.2. 新建栅格数据集¶

```driver = inDataset.GetDriver()
```

```Create(<filename>, <xsize>, <ysize>, [<bands>], [<GDALDataType>])
```

```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()
```

## 6.3. 计算band的统计量¶

```outBand.FlushCache()
outBand.GetStatistics(0, 1)
```

## 6.4. 设定新图的地理参考点¶

```geoTransform = inDataset.GetGeoTransform()
outDataset.SetGeoTransform(geoTransform )
proj = inDataset.GetProjection()
outDataset.SetProjection(proj)
```

## 6.5. 建立pyramids¶

```gdal.SetConfigOption('HFA_USE_RRD', 'YES')
```

```outDataset.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])
```

## 6.6. 图像的拼接¶

1. 对每张图：读取行数和列数，原点(minX,maxY)，像素长，像素宽，并计算坐标范围
```maxX1 = minX1 + (cols1 * pixelWidth)
minY1 = maxY1 + (rows1 * pixelHeight)
```
1. 计算输出图像的坐标范围：
```minX = min(minX1, minX2, …) maxX = max(maxX1, maxX2, …)
minY = min(minY1, minY2, …) maxY = max(maxY1, maxY2, …)
```
1. 计算输出图像的行数和列数：
```cols = int((maxX – minX) / pixelWidth)
rows = int((maxY – minY) / abs(pixelHeight)
```
1. 建立并初始化输出图像
2. 对每张待拼接的图：计算offset值
```xOffset1 = int((minX1 - minX) / pixelWidth)
yOffset1 = int((maxY1 - maxY) / pixelHeight)
```

6. 对输出图像：计算统计量，设定geotransform ：

```[minX, pixelWidth, 0, maxY, 0, pixelHeight]
```