4_获取影像像元数据

获取影像像元数据

本文中的方法可以对第三章进行补充，能够访问像元灰度值，并对其进行处理。

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)

  读取图像数据(以二进制的形式)

• 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)

  读取图像数据(以数组的形式)

参数说明：

xoff, yoff指定想要读取的部分原点位置在整张图像中距离全图原点的位置；

xsize, ysize指定要读取部分图像的矩形的长和宽；

buf_xsize, buf_ysize~~放缩比例，最后读取影像的长宽会根据这个比例因子进行缩放；~~测试发现这两个参数并不是缩放比例因子，暂时不清楚这两个参数的意义，但是官方文档中提到，通常不应该使用这两个参数。

---

已解决：

buf_xsize, buf_ysize并不是缩放比例因子，而是新的指定区域，比如指定xsize, ysize为(20, 20)，那么buf_xsize, buf_ysize设置为(50, 50)或者(10, 10)，就说明通过已经提取的20矩阵，计算成为50矩阵或者10矩阵。

因此：

进行分析运算时，只需要前4个参数即可，后面两个参数是用于将截取的影像强行放置在限定范围区域中进行缩放。

---

from osgeo import gdal
from osgeo.gdalconst import *

tif1 = gdal.Open('I:\\MasterLife\\资源与环境遥感\\实验作业一\\caijian2.tif')
print(tif1.ReadAsArray())  # 获取整幅影像的像元数组
help(tif1.ReadAsArray)  # 参看函数的帮助
print(tif1.ReadAsArray(xoff = 300, yoff = 300, xsize = 6, ysize = 6))
print(tif1.ReadRaster(xoff = 300, yoff = 300, xsize = 6, ysize = 6))

ReadRaster返回的结果是一个二进制的bytearray；

ReadAsArray返回结果是numpy的数组，类型为int16，方便与numpy结合进行分析和处理。

# ReadAsArray
[[[ 7  2  4  4  4  3]
  [ 3  5  5  5  6  3]
  [14 16 22 28 21 16]
  [40 51 42 35 32 38]
  [36 42 32 24 27 34]
  [28 24 23 19 24 17]]

 [[17  4  3  9  4  7]
  [ 5 10  9  7  9  5]
  [13 19 34 33 30 26]
  [60 67 60 56 49 53]
  [58 63 52 43 42 51]
  [45 42 44 35 44 27]]

 [[ 4  0  0  2  0  0]
  [ 0  2  0  6  4  0]
  [11 19 32 34 27 24]
  [59 67 55 48 44 51]
  [50 61 42 31 35 45]
  [34 29 30 23 32 16]]]

# ReadRaster
bytearray(b'\x07\x02\x04\x04\x04\x03\x03\x05\x05\x05\x06\x03\x0e\x10\x16\x1c\x15\x10(3*# &$* \x18\x1b"\x1c\x18\x17\x13\x18\x11\x11\x04\x03\t\x04\x07\x05\n\t\x07\t\x05\r\x13"!\x1e\x1a<C<815:?4+*3-*,#,\x1b\x04\x00\x00\x02\x00\x00\x00\x02\x00\x06\x04\x00\x0b\x13 "\x1b\x18;C70,32=*\x1f#-"\x1d\x1e\x17 \x10')

如果想要读取整个影像则：

# 读取整个影像
tif1.ReadAsArray()
tif1.ReadAsArray(xoff=0, yoff= 0,xsize=tif1.RasterXSize, ysize=tif1.RasterYSize)

注意：

Band和RasterSet下的ReadAsArray函数的参数有一定的区别。

区别在于：

在RasterSet在的两个参数是xsize, ysize，在Band中是win_xsize, win_ysize

获取单个波段

注意：这里GetRasterBand()是从1开始的，不是从0.

band1 = tif1.GetRasterBand(1)
print(band1.ReadAsArray())

获取范围越界

print(band1.XSize, band1.YSize)  # 2387 2637
print(band1.ReadAsArray(0, 0, 3000, 3000))
# ERROR 5: I:\MasterLife\资源与环境遥感\实验作业一\caijian2.tif
# band 1: Access window out of range in RasterIO().  Requested
# (0,0) of size 3000x3000 on raster of 2387x2637.

会进行ERROR 5的错误提示。

分块读取影像

测验分块横向、竖向读取和直接读取影像之间的速度差异。

dataset1 = gdal.Open('I:\\MasterLife\\资源与环境遥感\\实验作业一\\caijian2.tif')
dataset = gdal.Open('I:\\MasterLife\\遥感机理\结课作业\\N11_RPRO_AVH_L1B_1P_19920417T060234_19920417T060633_018353_v0100\\N11_RPRO_AVH_L1B_1P_19920417T060234_19920417T060633_018353\\N11_RPRO_AVH_L1B_1P_19920417T060234_19920417T060633_018353\\image.l1b')

band1 = dataset.GetRasterBand(1)
# band的长和宽
band1_x = band1.XSize
band1_y = band1.YSize
# 对影像分块
# 块宽
block_item = 128
# 需要先读取一遍，目的是进行缓存一次，否则第一次读取的速度会额外大
band1.ReadAsArray(0, 0, band1_x, band1_y)
# 第一次读取
start_time1 = time.time()
band1.ReadAsArray(0, 0, band1_x, band1_y)
start_time2 = time.time()
# 两个方向的块数
x_items = int(band1_x / block_item)
y_items = int(band1_y / block_item)
# 分块读取
# 纵向
start_time3 = time.time()
for i in range(x_items):
    for j in range(y_items):
        band1.ReadAsArray(128 * i, 128 * j, x_items, y_items)
start_time4 = time.time()
# 横向
for i in range(y_items):
    for j in range(x_items):
        band1.ReadAsArray(128 * j, 128 * i, x_items, y_items)
start_time5 = time.time()
print(f"直接读取：{start_time2 - start_time1}\n"
      f"纵向读取：{start_time4 - start_time3}\n"
      f"横向读取：{start_time5 - start_time4}")


直接读取：0.003996849060058594
纵向读取：0.002991199493408203
横向读取：0.0019958019256591797

地图代数

需要使用numpy包

有点小毛病，不知道怎么解决numpy的问题。

波段数据类型

即DataType，定义在gdalconst中，可以通过dir(gdalconst)查看该常数的所有属性和方法。基本上gdalconst中封装的全是常量。

其中GDT开头的就是全部的波段数据类型。

from osgeo import gdal, gdalconst
from osgeo.gdalconst import *

print(dir(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

访问索引图像的处理

这里的索引图像需要是有颜色表索引生成的影像，包含ColorMap、色彩空间等一系列概念，因暂时没有索引图像而跳过这部分。

2.5. 访问索引图像的处理 — 首页 (osgeo.cn)