【GDAL-Python】10-在Python中可视化多波段卫星影像

1-介绍

1.1 主要内容

（1）在本教程中，主要介绍如何使用 Python 和 matplotlib 可视化多波段 Landsat 8 卫星影像组成的真彩色影像以及假彩色合成影像。

（2）基本过程：
1）将单个栅格波段读取为数组；
2）显示图像时的三种显示策略将波段像素值缩放到0-1；
3）使用 numpy.dstack 堆叠数组 ；
4）使用 plt.imshow（） 绘制 RGB 合成

（3）视频地址：B站对应教程-在Python中可视化多波段卫星影像

（4）显示策略
由于卫星影像的像素值一般不为8bit（0-255），通常为16bit或者32bit；参考QGIS中显示图像的策略，本教程共三种显示策略，分别是：

• 1）归一化值显示：（像素值-最小值）/(最大值-最小值)
• 2）线性拉伸显示：（像素值-百分比累计值1）/（百分比累计值2-百分比累计值1）；
• 3）标准差值后显示
  PS：遥感图像中最常用的是线性拉伸。

参考QGIS显示策略如下图：

1.2 线性拉伸介绍

遥感卫星影像拉伸是一种图像处理技术，主要用于调整或增强图像的亮度和对比度，以更清晰地显示图像中的特征和细节，从而改善图像质量并帮助更准确地进行地理空间分析和解译，本质上是对像素值进行重新映射：具体技术见参考资料（1）。
2%的线性拉伸：选用像素值中2%的值都在百分比累计值阈值内的值为上述百分比累计值1，同理选用像素值中98%的值都在百分比累计值阈值内的值为上述百分比累计值2，以此两个值来计算重新映射后的像素值。

2-代码实现

2.1 数据介绍

（1）原始数据信息：选用了4个波段的WordView3影像，影像分辨率为1.24米，理论上有八个波段，但选取的实验数据为4个波段，影像尺寸为：影像元数据如下图：

（2）原始影像在QGIS显示效果：
​

2.2 代码实现

#！/usr/bin/env python3
# -*- coding:utf-8 -*-
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
"""1.归一化值显示函数"""
def scaleMinMax(x):
    return((x - np.nanmin(x))/(np.nanmax(x) - np.nanmin(x)))
"""2.2%线性拉伸显示函数"""
def scaleCCC(x):
    return((x - np.nanpercentile(x, 2))/(np.nanpercentile(x, 98) - np.nanpercentile(x,2)))
"""3.标准差值显示函数"""
def scaleStd(x):
    return((x - (np.nanmean(x)-np.nanstd(x)*2))/((np.nanmean(x)+np.nanstd(x)*2) - (np.nanmean(x)-np.nanstd(x)*2)))
#==========================（1）读取多光谱影像========================
ds = gdal.Open(r"GDAL_testing_data\JAX_IMG1_MSI.TIF")
#==========================（2）读取出对应rgb的对应波段================
r = ds.GetRasterBand(1).ReadAsArray()
g = ds.GetRasterBand(2).ReadAsArray()
b = ds.GetRasterBand(3).ReadAsArray()

ds = None
#==========================（3）归一化值处理并显示=====================
rMinMax = scaleMinMax(r)
gMinMax = scaleMinMax(g)
bMinMax = scaleMinMax(b)

rgbMinMax = np.dstack((rMinMax,gMinMax,bMinMax))
plt.figure()
plt.imshow(rgbMinMax)
plt.show()
#==========================（4）2%线性拉伸处理并显示===================
rCCC = scaleCCC(r)
gCCC = scaleCCC(g)
bCCC = scaleCCC(b)

rgbCCC = np.dstack((rCCC,gCCC,bCCC))
plt.figure()
plt.imshow(rgbCCC)
plt.show()

rStd = scaleStd(r)
gStd = scaleStd(g)
bStd = scaleStd(b)
#==========================（5）标准差值处理并显示=====================
rgbStd = np.dstack((rStd,gStd,bStd))
plt.figure()
plt.imshow(rgbStd)
plt.show()

2.3 效果显示

运行效果与QGIS打开显示效果一致。

4-参考资料

（1）遥感图像显示拉伸