Python遥感图像处理指南(3)-光谱分析

        这是Python遥感图像处理的的第三章，使用 Python 进行空间数据分析，并实现地理空间相关流程的自动化。在该系列的第一篇文章中，我们学习了如何使用 Anaconda 软件包管理器在 Windows 上准备 Python 环境，以及如何从 Jupyter 笔记本打开 GeoTiff 图像。接下来，在第二篇文章中，我们了解了矩阵操作的基础知识，还学习了如何创建一个灵活的函数，在给定两个波段的情况下计算归一化差异指数（如 NDVI、MNDWI 等）。

        在第三篇文章中，我们将继续分析图像数据，目标是绘制不同目标的光谱数据。今天的文章将指导利用Python显示水或植被的平均光谱曲线。

1、了解数组切片切片

        在第一篇文章中，我们看到了如何通过向数组变量传递方括号中的索引来访问单个像素值，例如：img[3000, 3000]。然而，访问单个值在卫星图像中并没有太大的用处，我们可以轻松处理数百万个值。因此，在进行光谱分析之前，我们必须学习的第一个概念是如何裁剪数组的一部分，这就是所谓的数组切片。NumPy 将 Python 的基本概念 "切片"（通常用于列表、字符串或元组）扩展到了 N 维数组，我们可以用它来访问图像的子集。它可以是一条线、一行、一个子数组或均匀分布的点。

        对数组进行切分的基本语法是在方括号中传递参数，参数之间用:分隔。例如，一个名为 v 的一维数组：v[first_index:last_index:step］，第一个元素的索引总是 0（这称为基于零的索引）；例如，如果我们有一个包含 5 个元素的一维数组，并希望抓取 3 个元素，则应使用 [0:3]，

        对于一个以上的维度（在我们的示例中，我们有 2 个或更多维度），我们应该为每个维度传递一个片段，片段之间用", "分隔。例如：img[0:1000, 0:1000]；索引可以倒数，例如使用-1（最后一个值）、-2、-3；当切片中的索引被省略时，表示我们需要该 "方向 "上的所有元素。例如：img[:2000, -3000:] 将访问前 2000 行（从 0 到 2000）和后 3000 行（从 -3000 到结束），如下图 所示。

2、放大显示某个区域

        下面的代码用于打开图像，城市大约位于第 4000 行和 5000 行像素之间，以及第 5400 行和 6400 行之间，因此我们将使用这些值进行切片。

import numpy as np
import matplotlib.pyplot as plt
import rasterio
from pathlib import Path

def load_landsat_image(img_folder, bands):
    image = {}
    path = Path(img_folder)
    for band in bands:
        # considering the landsat images end with *_SR_B#.TIF, we will use it to locate the correct file
        file = next(path.glob(f'*_SR_{band}.tif'))
        print(f'Opening file {file}')
        ds = rasterio.open(file)
        image.update({band: ds.read(1)})

    return image

# load the image
img = load_landsat_image('D:/Images/Landsat/LC08_L2SP_231062_20201026_20201106_02_T1/', ['B2', 'B3', 'B4'])

# stack the layers to create a cube
rgb = np.stack([img['B4'], img['B3'], img['B2']], axis=-1)

# normalize the values
rgb = rgb/rgb.max() * 2

# display the image with a slightly increased figure size
plt.f