>>> from env_helper import info; info()
页面更新时间: 2024-01-07 20:30:06
运行环境:
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-16-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

2.4. 使用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, buf_obj=None) method of osgeo.gdal.Dataset instance

上面的 help() 函数查看gdal方法。

继续查看 ReadAsArray() 方法。

>>> 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, interleave='band', band_list=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

这是两个非常重要的函数,它们可以直接读取图像的数据,从而对栅格数据进行分析。

可以看到两个函数有许多的参数。解释一下:

  • xoff,yoff :指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位)。

  • xsize,ysize : 指定要读取部分图像的矩形的长和宽(以像元为单位)。

  • buf_xsize,buf_ysize :可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高, GDAL 将帮你缩放到这个大小。

  • buf_type :可以对读出的数据的类型进行转换(比如原图数据类型是short,你要把它们缩小成byte)。

  • band_list :适应多波段的情况。可以指定要读取的波段。

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

>>> 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)
bytearray(b'x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00x0cx00')

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

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

2.4.1. 读取波段中的数据

ReadRaster()的函数细节

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

>>> from osgeo import gdal
>>> # from osgeo.gdalconst import *
>>> dataset = gdal.Open("/gdata/geotiff_file.tif")
>>> band = dataset.GetRasterBand(1)
>>> band.ReadAsArray(100,100,5,5,10,10)
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参数的作用就是把硬盘上指定区域 (前四个参数定义)的数据按比例缩放到用户指定区域(后两个定义)内, 必要的时候进行缩放。 这里需要注意的是缩的情况(缩的时候是取几个周围点的平均值) 如果是调色板数据就可能引发问题。

栅格数据范围的处理

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

>>> band.XSize
1500
>>> band.YSize
900
>>> band.ReadAsArray(95,100,5,5)
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模块可以很方便地扩展矩阵。

2.4.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.018355846405029297
0.059012413024902344

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

运行一下得到结果:

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

0.03625917434692383
0.06595683097839355
0.06917715072631836

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

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

2.4.3. 地图代数

以计算NDVI为例:

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

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

编程要点如下:

  1. 将波段3读入数组 data3 ,将波段2读入数组 data2

  2. 使用上面的计算公式;

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

代码如下:

>>> 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)))
/tmp/ipykernel_62017/3808589461.py:11: RuntimeWarning: invalid value encountered in 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)

2.4.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与整型的对应值。

  • 未知或未指定类型 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