GDAL python教程基础篇（10）GDAL读取数据集数据

访问遥感影像的描述性信息，可以概括地知道影像的获取时间、处理时间、空间分辨率、影像大小等一些信息。 但是为了对遥感影像进行处理，需要进一步访问遥感影像中的数据，即影像中像元的灰度值。

GDAL提供了下面两个函数来访问影像的数值。
ReadRaster() 读取图像数据(以二进制的形式)
ReadAsArray() 读取图像数据(以数组的形式)

>>> from osgeo import gdal
>>> dataset = gdal.Open("/gdata/lu75c.tif")
>>> help(dataset.ReadRaster)
>Help on method ReadRaster in module osgeo.gdal:

ReadRaster(xoff=0, yoff=0, xsize=None, ysize=None, buf_xsize=None, buf_ysize=None, buf_type=None, band_list=None, buf_pixel_space=None, buf_line_space=None, buf_band_space=None, resample_alg=0, callback=None, callback_data=None) method of osgeo.gdal.Dataset instance

上面的 help() 函数查看gdal方法:退出查看页面在终端输入“q”。
这是两个非常重要的函数，它们可以直接读取图像的数据， 从而对栅格数据进行分析。可以看到两个函数的帮助中有许多的参数。 解释一下：
xoff,yoff ：指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位）。
xsize,ysize ： 指定要读取部分图像的矩形的长和宽（以像元为单位）。
buf_xsize，buf_ysize ：可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高，gdal将帮你缩放到这个大小。
buf_type ：可以对读出的数据的类型进行转换（比如原图数据类型是short，你要把它们缩小成byte）。
band_list ：适应多波段的情况。可以指定要读取的波段。

>>> import array
>>> from numpy import *
>>> dataset.ReadAsArray(2500,2500,3,3)
>>> array([[12, 12, 12],
>>> [12, 12, 12],
>>> [12, 12, 12]],dtype=int16)
>>> dataset.ReadRaster(2500,2500,3,3)
b'x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00'

ReadAsArray()读出的是numpy的数组，类型为int16； 而ReadData()读出的是二进制型。

1.读取波段中的数据

ReadAsArray() ，返回的是numpy模块定义的Array， 之所以用它的是因为它排列很工整。

from osgeo import gdal
from osgeo.gdalconst import *
dataset = gdal.Open(r"D:\work\遥感图像处理\数据\BOA Reflectance-10m_MTD_MSIL2A.tif")
band = dataset.GetRasterBand(1)
print(band.ReadAsArray(100,100,5,5,10,10))

前两个100是取值窗口的左上角在实际数据中所处象元的(x, y)位置。 后两个是取值窗口覆盖的区域大小， 再后面是取值窗口取出数组进行缩放后数组的大小。 这里需要注意的是这里的buffer大小是根据参数自动分配的，可以不指定， 如果不指定，则和第3，4个参数一致。经过5，6两个参数的设置，可以进行缩放。 假如取值窗口大小是20 × 20，取出后数组就可以人工设置大小。 让它成为10 × 10的数组，或者是40 × 40的数组。 如果设置成20 × 40的数组则取出的数组对于真实图像来说有了变形。

2.读取栅格数据方式与效率

对于GeoTiff来说， 从横向读取和纵向读取的效率相差很大。 可以写一个Python脚本文件来验证一下：

>>> from osgeo import gdal
>>> import time
>>> dataset = gdal.Open("/gdata/lu75c.tif")
>>> band = dataset.GetRasterBand(1)
>>> width, height = dataset.RasterXSize, dataset.RasterYSize
>>> bw, bh = 128, 128
>>> bxsize = width/bw
>>> bysize = height/bh
>>> band.ReadAsArray(0,0,width,height)
>>> start = time.time()
>>> band.ReadAsArray(0,0,width,height)
>>> print (time.time()-start)
>>> start2 = time.time()
>>> for i in range(int(bysize)):
>>>     for j in range(int(bxsize)):
>>>         band.ReadAsArray(bw*j,bh*i,bw,bh)
>>> print (time.time()-start2)
0.030247211456298828
0.06247711181640625

3. 地图代数

以计算NDVI为例：

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

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

编程要点如下：
将波段3读入数组 data3 ，将波段2读入数组 data2 ；
使用上面的计算公式；
当 data3 和 data2 均为 0 时（例如用 0 表示 NODATA ）， 会出现被0除的错误，导致程序崩溃。需要用 mask 配合 choose 将 0 值去掉。

代码如下：

>>> from osgeo import gdal
>>> import numpy
>>> dataset = gdal.Open("/gdata/geotiff_file.tif")
>>> band2 = dataset.GetRasterBand(2)
>>> band3 = dataset.GetRasterBand(3)
>>> cols = 100
>>> rows = 100
>>> data2 = band2.ReadAsArray(0, 0, cols,rows).astype(numpy.float16)
>>> data3 = band3.ReadAsArray(0, 0, cols,rows).astype(numpy.float16)
>>> mask = numpy.greater(data3 + data2, 0)
>>> ndvi = numpy.choose(mask, (-99, (data3 - data2) / (data3 + data2)))
<ipython-input-10-2e95808c0843>:11: RuntimeWarning: invalid value encountered in true_divide
  ndvi = numpy.choose(mask, (-99, (data3 - data2) / (data3 + data2)))
>>> ndvi
array([[0.2766, 0.318 , 0.4285, ..., 0.3057, 0.275 , 0.6   ],
       [0.1951, 0.2208, 0.3135, ..., 0.2683, 0.3538, 1.    ],
       [0.2142, 0.2195, 0.2896, ..., 0.4182, 0.758 , 1.    ],
       ...,
       [0.6313, 0.6113, 0.6553, ..., 0.3333, 0.359 , 0.36  ],
       [0.4119, 0.4119, 0.422 , ..., 0.3784, 0.3784, 0.3418],
       [0.3684, 0.322 , 0.309 , ..., 0.4084, 0.4   , 0.3684]],
      dtype=float16)

4. 波段数据

DataType是指图像中实际数值的数据类型。具体数据类型定义在gdalconst模块里。
使用的时候用from osgeo import gdalconst引入。
返回结果如下：

>>> from osgeo import gdalconst
>>> dir(gdalconst)[:6] + ['... ...'] + dir(gdalconst)[-3:]
['CE_Debug',
 'CE_Failure',
 'CE_Fatal',
 'CE_None',
 'CE_Warning',
 'CPLES_BackslashQuotable',
 '... ...',
 '_swig_repr',
 '_swig_setattr_nondynamic_class_variable',
 '_swig_setattr_nondynamic_instance_variable']

那些GDT开头的就是数值数据类型。

想要查看图像中某一波段的数据类型，只需要打印这一波段的DataType属性即可。

>>> from osgeo import gdal
>>> dataset = gdal.Open("/gdata/geotiff_file.tif")
>>> band = dataset.GetRasterBand(1)
>>> band.DataType
1

返回结果为整型。原来1表示的是 gdalconst.GDT_Byte 。注意这里的类型是与numpy中的类型对应的。下面我们来看一个gdalconst与整型的对应值。