文章目录
- 1-介绍
- 1.1 主要内容
- 1.2 线性拉伸介绍
- 2-代码实现
- 2.1 数据介绍
- 2.2 代码实现
- 2.3 效果显示
- 4-参考资料
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)遥感图像显示拉伸