python矩阵堆叠-实现遥感影像波段叠加

前文读取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)