GDAL python教程（5）——地图代数与栅格数据的写入
====================================================
    

以计算NDVI为例：
-------------------------------------------------
::

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

其中NIR为波段3，RED为波段2

编程要点如下：

1. 将波段3读入数组data3，将波段2读入数组data2
2. 计算公式为：  ndvi = (data3 – data2) / (data3 + data2)
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(<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()

 

计算band的统计量
-----------------------------------------

首先用FlushCache()把缓存数据写入磁盘

之后用GetStatistics(<approx_ok>, <force>)计算统计量。如果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