6.1. 以计算NDVI为例:
NDVI=(NIR-RED)/(NIR+RED)
其中NIR为波段3,RED为波段2
编程要点如下:
- 将波段3读入数组data3,将波段2读入数组data2
- 计算公式为: ndvi = (data3 – data2) / (data3 + data2)
- 当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)))
6.2. 新建栅格数据集
将刚才计算得到的数据写入新的栅格数据集之中
首先要复制一份数据驱动:
driver = inDataset.GetDriver()
之后新建数据集
Create(<filename>, <xsize>, <ysize>, [<bands>], [<GDALDataType>])
其中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()
6.3. 计算band的统计量
首先用FlushCache()把缓存数据写入磁盘
之后用GetStatistics(<approx_ok>, <force>)计算统计量。如果approx_ok=1那么计算是基于pyramid的,如果force=0那么当整幅图都要被重读一遍的时候就不计算统计量了。
outBand.FlushCache() outBand.GetStatistics(0, 1)
6.4. 设定新图的地理参考点
如果新图与另一张图的地理参考信息完全一致,那就很简单了
geoTransform = inDataset.GetGeoTransform() outDataset.SetGeoTransform(geoTransform ) proj = inDataset.GetProjection() outDataset.SetProjection(proj)
6.5. 建立pyramids
设定Imagine风格的pyramids
gdal.SetConfigOption('HFA_USE_RRD', 'YES')
强制建立pyramids
outDataset.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])
6.6. 图像的拼接
- 对每张图:读取行数和列数,原点(minX,maxY),像素长,像素宽,并计算坐标范围
maxX1 = minX1 + (cols1 * pixelWidth) minY1 = maxY1 + (rows1 * pixelHeight)
- 计算输出图像的坐标范围:
minX = min(minX1, minX2, …) maxX = max(maxX1, maxX2, …) minY = min(minY1, minY2, …) maxY = max(maxY1, maxY2, …)
- 计算输出图像的行数和列数:
cols = int((maxX – minX) / pixelWidth) rows = int((maxY – minY) / abs(pixelHeight)
- 建立并初始化输出图像
- 对每张待拼接的图:计算offset值
xOffset1 = int((minX1 - minX) / pixelWidth) yOffset1 = int((maxY1 - maxY) / pixelHeight)
读入数据并按照上面计算的offset写入
6. 对输出图像:计算统计量,设定geotransform :
[minX, pixelWidth, 0, maxY, 0, pixelHeight]
设定projection,建立pyramids