基于Python 进行卫星图像多种指数分析

一、前言

本文帮助读者更好地了解卫星数据以及使用 Python 探索和分析哨兵2卫星数号数据在Sundarbans地区的不同方法。

二、Sundarbans研究区

孙德尔本斯（Sundarbans）是恒河、雅鲁藏布江和梅克纳河在孟加拉湾汇合形成的三角洲中最大的红树林区之一。 孙德尔本斯森林横跨印度和孟加拉国，面积约 10,000 平方公里，其中 40% 位于印度境内，是许多珍稀和全球受威胁野生动物物种的家园。 下面的谷歌地图显示了孙德尔本斯地区。

在本文中，我们将使用2020年1月27日使用Sentinel-2卫星获取的Sundarbans卫星数据的一小部分。该Sundarbans数据的优化自然色图像如下所示：

卫星数据有 954 * 298 像素，12 个波段，光谱分辨率从 10 — 60 米不等。

三、数据分析

为了对 Sundarbans 数据执行不同的操作，我们将使用 EarthPy、RasterIO、Matplotlib、Plotly 等库进行数据可视化和分析。

让我们通过导入必要的包开始编码：

from glob import globimport earthpy as et
import earthpy.spatial as es
import earthpy.plot as epimport rasterio as rioimport matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import ListedColormapimport plotly.graph_objects as gonp.seterr(divide='ignore', invalid='ignore')

（1）数据读取

让我们使用 rasterio 读取 12 个波段，并使用 numpy.stack() 方法将它们堆叠到一个 n 维数组中。 堆叠后的结果数据具有形状 (12, 954, 298)。

S_sentinel_bands = glob("/content/drive/MyDrive/Satellite_data/sundarbans_data/*B?*.tiff")S_sentinel_bands.sort()l = []for i in S_sentinel_bands:with rio.open(i, 'r') as f:l.append(f.read(1))arr_st = np.stack(l)

（2）波段可视化

正如我们所讨论的，数据包含 12 个波段。 让我们使用 EarhPy 包可视化每个波段。 plot_bands() 方法获取波段和绘图的堆栈以及自定义标题，这可以通过使用 title= 参数将每个图像的唯一标题作为标题列表传递来完成。

ep.plot_bands(arr_st, cmap = 'gist_earth', figsize = (20, 12), cols = 6, cbar = False)
plt.show()

（3）RGB波段合成

研究区有多个波段，包含从可见光到红外线的数据。 因此很难为人类可视化数据。 通过创建 RGB 合成图像，可以更轻松有效地理解数据。 要绘制 RGB 合成图像，您将绘制红色、绿色和蓝色波段，它们分别是波段 4、3 和 2。 由于 Python 使用从零开始的索引系统，因此您需要从每个索引中减去 1 的值。 因此，红色波段的索引为 3，绿色为 2，蓝色为 1。

如果像素亮度值偏向零值，我们创建的合成图像有时会变暗。 此类问题可以通过拉伸图像中的像素亮度值来解决，使用参数 stretch=True 将值扩展到潜在值的完整 0-255 范围，以增加图像的视觉对比度。 此外，str_clip 参数允许您指定要剪掉多少数据尾部。 数字越大，数据被拉伸或变亮的程度就越大。

让我们看看绘制 RGB 合成图像以及应用拉伸的代码。

# RGB Composite Imagergb = ep.plot_rgb(arr_st, rgb=(3,2,1), figsize=(10, 16))
plt.show()# RGB Composite Image with Strechep.plot_rgb(arr_st,rgb=(3, 2, 1),stretch=True,str_clip=0.2,figsize=(10, 16))
plt.show()

上图显示了应用拉伸前后孙德尔本斯数据的 RGB 合成图像。

（4）直方图

可视化高光谱图像数据集的波段有助于我们了解波段像素/值的分布。 earhtpy.plot 中的 hist 方法通过为我们之前创建的数据集/堆栈的波段绘制直方图来完成这项工作。 我们还可以修改各个直方图的列大小、标题和颜色。 让我们看看绘制直方图的代码。

colors = ['tomato', 'navy', 'MediumSpringGreen', 'lightblue', 'orange', 'blue','maroon', 'purple', 'yellow', 'olive', 'brown', 'cyan']ep.hist(arr_st, colors = colors,title=[f'Band-{i}' for i in range(1, 13)], cols=3, alpha=0.5, figsize = (12, 10))plt.show()

（5）植被和土壤指数

归一化卫星指数是根据多光谱卫星图像计算得出的图像。 这些图像强调存在的特定现象，同时减轻降低图像效果的其他因素。 例如，植被指数将在索引图像中将健康的植被显示为明亮，而不健康的植被具有较低的值，贫瘠的地形是黑暗的。 由于地形变化（山丘和山谷）的阴影会影响图像的强度，因此创建指数的方式是强调对象的颜色而不是对象的强度或亮度。

• 归一化植被指数 (NDVI)

要确定一块土地上的绿色密度，研究人员必须观察植物反射的可见光 (VIS) 和近红外 (NIR) 阳光的不同颜色（波长）。 归一化差异植被指数 (NDVI) 通过测量植被强烈反射的近红外光和红光（植被吸收的）之间的差异来量化植被。 NDVI 的范围总是从 -1 到 +1。

NDVI = ((NIR - Red)/(NIR + Red))

例如，当你有负值时，它很可能是水。 另一方面，如果您的 NDVI 值接近 +1，则它很可能是茂密的绿叶。 但当 NDVI 接近于零时，就没有绿叶，甚至可能是城市化区域。

让我们看看在 Sundarbans 卫星数据上实现 NDVI 的代码。

ndvi = es.normalized_diff(arr_st[7], arr_st[3])ep.plot_bands(ndvi, cmap="RdYlGn", cols=1, vmin=-1, vmax=1, figsize=(10, 14))plt.show()

• 土壤调整植被指数 (SAVI)

土壤调整植被指数 (SAVI) 是一种植被指数，它试图使用土壤亮度校正因子来最小化土壤亮度的影响。 这通常用于植被覆盖率低的干旱地区。

SAVI = ((NIR - Red) / (NIR + Red + L)) x (1 + L)

L 值根据绿色植被覆盖量而变化。 一般来说，在没有绿色植被覆盖的地区，L=1； 在中度绿色植被覆盖区，L=0.5； 在植被覆盖度很高的地区，L=0（相当于NDVI方法）。 该索引输出介于 -1.0 和 1.0 之间的值。 让我们看看SAVI的实现代码。

L = 0.5savi = ((arr_st[7] - arr_st[3]) / (arr_st[7] + arr_st[3] + L)) * (1 + L)ep.plot_bands(savi, cmap="RdYlGn", cols=1, vmin=-1, vmax=1, figsize=(10, 14))plt.show()

• 可见大气阻力指数 (VARI)

可见大气阻力指数 (VARI) 旨在强调光谱可见部分的植被，同时减轻光照差异和大气影响。 它是 RGB 或彩色图像的理想选择； 它利用了所有三个色带。

VARI = (Green - Red)/ (Green + Red - Blue)

vari = (arr_st[2] - arr_st[3])/ (arr_st[2] + arr_st[3] - arr_st[1])ep.plot_bands(vari, cmap="RdYlGn", cols=1, vmin=-1, vmax=1, figsize=(10, 14))plt.show()

（6）水体指数

地表水变化是环境、气候和人类活动的一个非常重要的指标。 遥感器，例如 sentinel-2、Landsat，在过去四十年中一直在提供数据，这些数据对于提取土地覆盖类型（例如森林和水）很有用。 研究人员提出了许多地表水提取技术，其中基于指数的方法因其简单性和成本效益而广受欢迎。

• 修正归一化差水指数 (MNDWI)

修正归一化差异水域指数 (MNDWI) 使用绿色和 SWIR 波段来增强开放水域特征。 它还减少了在其他指数中通常与开阔水域相关的建成区特征。

MNDWI = (Green - SWIR) / (Green + SWIR)

下面的代码用于实现 MNDWI，输出如下所示。

mndwi = es.normalized_diff(arr_st[2], arr_st[10])ep.plot_bands(mndwi, cmap="RdYlGn", cols=1, vmin=-1, vmax=1, figsize=(10, 14))plt.show()

• 归一化差异水分指数 (NDMI)

归一化差异水分指数 (NDMI) 对植被中的水分含量很敏感。 它用于监测干旱以及监测火灾多发地区的燃料水平。 它使用 NIR 和 SWIR 波段来创建旨在减轻照明和大气影响的比率。

NDMI = (NIR - SWIR1)/(NIR + SWIR1)

让我们看看实现和输出：

ndmi = es.normalized_diff(arr_st[7], arr_st[10])ep.plot_bands(ndmi, cmap="RdYlGn", cols=1, vmin=-1, vmax=1, figsize=(10, 14))plt.show()

（6）地质指数

事实证明，卫星图像和航空摄影是支持矿产勘探项目的重要工具。 它们可以以多种方式使用。 首先，它们为地质学家和现场工作人员提供轨道、道路、围栏和居住区的位置。

• 粘土矿物

Clay Minerals Ratio = SWIR1 / SWIR2

cmr = np.divide(arr_st[10], arr_st[11])ep.plot_bands(cmr, cmap="RdYlGn", cols=1, vmin=-1, vmax=1, figsize=(10, 14))plt.show()

• 黑色金属矿产

黑色金属矿物比例突出了含铁材料。 它使用 SWIR 波段和 NIR 波段之间的比率。

Ferrous Minerals Ratio = SWIR / NIR

fmr = np.divide(arr_st[10], arr_st[7])ep.plot_bands(fmr, cmap="RdYlGn", cols=1, vmin=-1, vmax=1, figsize=(10, 14))plt.show()

四、结论

本文介绍了数据可视化和归一化植被、水和地质指数等不同方法，以使用 python 分析 Sundarbans 研究区卫星数据。