访问数据集的数据

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

GDAL提供了下面两个函数来访问影像的数值。

  • ReadRaster() 读取图像数据(以二进制的形式)
  • ReadAsArray() 读取图像数据(以数组的形式)
In [1]:
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”。继续查看 ReadAsArray() 方法。

In [2]:
help(dataset.ReadAsArray)
Help on method ReadAsArray in module osgeo.gdal:

ReadAsArray(xoff=0, yoff=0, xsize=None, ysize=None, buf_obj=None, buf_xsize=None, buf_ysize=None, buf_type=None, resample_alg=0, callback=None, callback_data=None) method of osgeo.gdal.Dataset instance
    Reading a chunk of a GDAL band into a numpy array. The optional (buf_xsize,buf_ysize,buf_type)
    parameters should generally not be specified if buf_obj is specified. The array is returned

同样需要键入“q”退出。

这是两个非常重要的函数,它们可以直接读取图像的数据, 从而对栅格数据进行分析。可以看到两个函数的帮助中有许多的参数。 解释一下:

  • xoff,yoff :指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位)。
  • xsize,ysize : 指定要读取部分图像的矩形的长和宽(以像元为单位)。
  • buf_xsize,buf_ysize :可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高,gdal将帮你缩放到这个大小。
  • buf_type :可以对读出的数据的类型进行转换(比如原图数据类型是short,你要把它们缩小成byte)。
  • band_list :适应多波段的情况。可以指定要读取的波段。

这里简单看一下如何获取GeoTIFF文件中的数据。

In [3]:
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)
Out[3]:
b'\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00'

我们就把图像左上角位于(3340, 3780),宽高都为10个像元的数据读取出来了。 这两个函数返回的结果不一样,其中ReadAsArray()读出的是numpy的数组,类型为int16; 而ReadData()读出的是二进制型,具体的解释见后面。

ReadAsArray()返回的结果。关于类型的更多解释,见下一节的表[tab:gdalconst]。

读取波段中的数据

ReadRaster()的函数细节

下面是我的一个例子:使用了 ReadAsArray() ,返回的是numpy模块定义的Array, 之所以用它的是因为它排列很工整。 (这里的数据是用GRASS从spearfish示例数据转出的aspect数据集)

In [4]:
from osgeo import gdal
from gdalconst import *
dataset = gdal.Open("/gdata/geotiff_file.tif")
band = dataset.GetRasterBand(1)
band.ReadAsArray(100,100,5,5,10,10)  
Out[4]:
array([[236, 236, 237, 237, 237, 237, 237, 237, 227, 227],
       [236, 236, 237, 237, 237, 237, 237, 237, 227, 227],
       [235, 235, 232, 232, 233, 233, 234, 234, 225, 225],
       [235, 235, 232, 232, 233, 233, 234, 234, 225, 225],
       [242, 242, 235, 235, 232, 232, 233, 233, 224, 224],
       [242, 242, 235, 235, 232, 232, 233, 233, 224, 224],
       [254, 254, 244, 244, 238, 238, 237, 237, 229, 229],
       [254, 254, 244, 244, 238, 238, 237, 237, 229, 229],
       [246, 246, 246, 246, 248, 248, 250, 250, 235, 235],
       [246, 246, 246, 246, 248, 248, 250, 250, 235, 235]], dtype=uint8)

ReadAsArray的原型

help(band.ReadAsArray)

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

band.ReadAsArray(100,100,10,10)

通过几个例子可以看到,读取的4个size参数的作用就是把硬盘上指定区域 (前四个参数定义)的数据按比例缩放到用户指定区域(后两个定义)内, 必要的时候进行缩放。 这里需要注意的是缩的情况(缩的时候是取几个周围点的平均值) 如果是调色板数据就可能引发问题。

栅格数据范围的处理

现在还有一个问题。就是取值窗口超过实际数据最大的范围怎么办?

In [5]:
band.XSize
Out[5]:
1500
In [6]:
band.YSize
Out[6]:
900
In [7]:
band.ReadAsArray(95,100,5,5)   
Out[7]:
array([[246, 242, 237, 233, 234],
       [247, 242, 237, 234, 234],
       [244, 244, 240, 239, 242],
       [241, 249, 245, 248, 254],
       [244, 251, 252, 251, 247]], dtype=uint8)
ERROR 5: Access window out of range in RasterIO().  Requested
(95,100) of size 5x5 on raster of 100x100.

可以看到,出错了。 获取数据的时候不能越界, 调用的时候要判断越界没。 还好用Python的numpy模块可以很方便地扩展矩阵。

读取栅格数据方式与效率

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

In [8]:
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.04159808158874512
0.06883883476257324

最后一段的循环的两个for位置调换一下。

运行一下得到结果:

运行python ch03_test_read_tif_time.py,得到以下结果:

0.03625917434692383
0.06595683097839355
0.06917715072631836

上面第9行,目的是为了先读取一遍数据,不然,第一次的结果会略大一些,可能会超过第二次,会让人以为把图像分块读取比不分块读取的效率要高。

至于出现这种情况的原因,在于影像的存储是有不同方式的。

地图代数

以计算NDVI为例:

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

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

编程要点如下:

  1. 将波段3读入数组data3,将波段2读入数组data2
  2. 计算公式为:

3.当data3和data2均为0时(例如用0表示NODATA),会出现被0除的错误,导致程序崩溃。需要用mask配合choose将0值去掉。

代码如下:

In [9]:
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)))
/home/bk/.local/lib/python3.5/site-packages/ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in true_divide
  # This is added back by InteractiveShellApp.init_path()
In [10]:
ndvi
Out[10]:
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)

波段数据类型

DataType是指图像中实际数值的数据类型。具体数据类型定义在gdalconst模块里。

使用的时候用from osgeo import gdalconst引入。

返回结果如下:

In [11]:
from osgeo import gdalconst
dir(gdalconst)[:6] + ['... ...'] + dir(gdalconst)[-3:]
Out[11]:
['CE_Debug',
 'CE_Failure',
 'CE_Fatal',
 'CE_None',
 'CE_Warning',
 'CPLES_BackslashQuotable',
 '... ...',
 '_swig_repr',
 '_swig_setattr',
 '_swig_setattr_nondynamic']

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

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

In [12]:
from osgeo import gdal
dataset = gdal.Open("/gdata/geotiff_file.tif")
band = dataset.GetRasterBand(1)
band.DataType
Out[12]:
1

返回结果为整型。原来1表示的是 gdalconst.GDT_Byte 。注意这里的类型是与numpy中的类型对应的。

下面我们来看一个gdalconst与整型的对应值。

  • 未知或未指定类型 gdalconst.GDT_Unknown 0
  • 8位无符整型 gdalconst.GDT_Byte 1
  • 16位无符整型 gdalconst.GDT_UInt16 2
  • 16位整型 gdalconst.GDT_Int16 3
  • 32位无符整型 gdalconst.GDT_UInt32 4
  • 32位整型值 gdalconst.GDT_Int32 5
  • 32位浮点型 gdalconst.GDT_Float32 6
  • 64位浮点型 gdalconst.GDT_Float64 7
  • 16位复数整型 gdalconst.GDT_CInt16 8
  • 32位复数整型 gdalconst.GDT_CInt32 9
  • 32位复数浮点型 gdalconst.GDT_CFloat32 10
  • 64位复数浮点型 gdalconst.GDT_CFloat64 11