前文读取landsat8数据,计算NDVI并保存成tiff文件和python批量读取landsat8的波段分别介绍了对landsat8数据的不同读取方式,以及对利用波段数据进行简单的指数计算。本文在此基础上实现波段的叠加。
准备工作
这一部分参考python批量读取landsat8的波段,利用文中的函数,我们对landsat8的波段路径进行了一个整理,存放到了bands列表中并返回。
看一下格式:第一张影像的前5个波段
from read_landsat8 import read_landsat8_bands
import numpy as np
from osgeo import gdal_array
base_path = 'data'
bands = read_landsat8_bands(base_path)
bands[0][:5]
[‘data\LC08_L1TP_121034_20130402_20200912_02_T1\LC08_L1TP_121034_20130402_20200912_02_T1_B1.TIF’, ‘data\LC08_L1TP_121034_20130402_20200912_02_T1\LC08_L1TP_121034_20130402_20200912_02_T1_B2.TIF’, ‘data\LC08_L1TP_121034_20130402_20200912_02_T1\LC08_L1TP_121034_20130402_20200912_02_T1_B3.TIF’, ‘data\LC08_L1TP_121034_20130402_20200912_02_T1\LC08_L1TP_121034_20130402_20200912_02_T1_B4.TIF’, ‘data\LC08_L1TP_121034_20130402_20200912_02_T1\LC08_L1TP_121034_20130402_20200912_02_T1_B5.TIF’]
读取波段1、2、3:
# 读取波段
B1 = bands[0][0]
B2 = bands[0][1]
B3 = bands[0][2]
B1_gdal = gdal_array.LoadFile(B1)
B2_gdal = gdal_array.LoadFile(B2)
B3_gdal = gdal_array.LoadFile(B3)
print(B1_gdal.shape)
print(B2_gdal.shape)
print(B3_gdal.shape) # 维度相等
(7301, 7341)
(7301, 7341)
(7301, 7341)
波段叠加:
本文只是简单演示,这里我们只对前面读取的123波段进行叠加
先用numpy转化成ndarray形式
# 转化成ndarray形式
B1_np = np.array(B1_gdal)
B2_np = np.array(B2_gdal)
B3_np = np.array(B3_gdal)
print(B1_np.shape)
print(B2_np.shape)
print(B3_np.shape)
维度仍保持
(7301, 7341)
(7301, 7341)
(7301, 7341)
用构建数组的方式叠加
B123 = np.array([B1_np, B2_np, B3_np])
print(B123.shape)
输出可以看成是channel_first形式
(3, 7301, 7341)
用stack方式叠加
B123 = np.stack([B1_np, B2_np, B3_np])
print(B123.shape)
输出维度与上面相同
(3, 7301, 7341)
stack叠加可以按照不同的轴实现不同的维度叠加,默认axis=0
B123 = np.stack([B1_np, B2_np, B3_np], axis = 0)
print(B123.shape)
(3, 7301, 7341)
B123 = np.stack([B1_np, B2_np, B3_np], axis = 1)
print(B123.shape)
(7301, 3, 7341)
B123 = np.stack([B1_np, B2_np, B3_np], axis = 2)
print(B123.shape)
(7301, 7341,3)
用concatenate方式叠加
B123 = np.concatenate([B1_np, B2_np, B3_np])
print(B123.shape)
可以看到,它是按照行叠加了
(21903, 7341)
这种叠加方式也有一些用处,例如,如果我们要认为构建一个mask矩阵用于存放标签。
# 获取宽、高
height = B123.shape[0]
width = B123.shape[1]
# 构建0-1 mask矩阵
mask = np.random.randint(0, 2, (height, width, 1))
print(mask.shape)
# 堆叠
img = np.concatenate([B123, mask], axis = 2)
print(img.shape)
(7301, 7341, 1)
(7301, 7341, 4)