python - 使用 python 库 Rasterio 创建 TIFF 图像的子集，然后显示并保存？

Question

我有两张光栅图像，一张来自波段 4，最后是 B4，另一张来自波段 5，最后是 B5。我想将 B5 栅格子集化为 800x600，然后显示它并将其保存为 GeoTiff。然后我想计算 NDVI（我假设我需要 B4 和 B5 来执行此操作，但不确定）。然后我想显示 B5 栅格的 NDVI 子集。显示它并将其另存为 GeoTiff。

我将如何创建 TIFF 光栅图像的 800 x 600 像素子集？我还想获取该 TIFF 并为该子集生成 NDVI 图像。

注意：我正在使用 Landsat 图像。图像在文件标题的末尾有 B5。

到目前为止我所做的：

import rasterio
from rasterio.windows import Window 
import matplotlib.pyplot plt # for later use

with rasterio.open('MyRasterImage.tif') as src:
    w = src.read(1, window=Window(0, 0, 800, 600))

我想使用 Spyder 或 Jupyter 笔记本来显示它。所以我想使用 matplotlib 并做了流动的代码：

# Plot
plt.imshow(w)
plt.show()

这样做会生成一个 800x600 matplotlib 窗口，但它都是紫色的，不知道为什么它会产生这个。

现在我希望能够显示这个 800x600 的图像。然后，我想在该子集 800x600 图像上执行 NDVI。然后显示带有 NDVI 显示的子集 800x600 图像。

我知道公式是：NDVI = (NIR - red) / (NIR + red)

但是我如何从这张 Landsat 图像中提取 NIR 和红色呢？

我的尝试：

band1 = dataset.read(1)
band2 = dataset.read(2)
band3 = dataset.read(3)
print(band[2])

当我为乐队运行该代码时，我收到错误：

rasterio indexerror: band index 2 out of range (not in (1,))

当我运行此代码时：

print(w.count)

它返回“1”。

所以这意味着陆地卫星图像只有一个波段？但是为了做 NDVI，我不需要 3 个波段吗？

我正在考虑编写一些这样的代码来从该栅格中获取 NDVI。但不知道如何去提取乐队：

# We handle the connections with "with"
with rasterio.open(bands[0]) as src:
    b3 = src.read(1)

with rasterio.open(bands[1]) as src:
    b4 = src.read(1)

# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

# Calculate NDVI
ndvi = (b4.astype(float) - b3.astype(float)) / (b4 + b3)

这段代码不起作用，因为波段没有被定义为任何东西，所以我不知道如何定义波段来获取 NDVI。

在此之后，我不确定如何显示和保存渲染图像。

标签： pythonrasterio