¿Crear un subconjunto de imágenes TIFF y generar su NDVI en Python?

Question

Tengo dos imágenes Raster, una de banda4 con un B4 al final y otra de banda5 con B5 al final. Quiero subconjuntar el raster B5 a 800x600, luego visualizarlo y guardarlo como un GeoTiff. Luego quiero calcular el NDVI (supongo que necesitaré tanto el B4 como el B5 para hacerlo, pero no estoy seguro). A continuación, quiero mostrar el subconjunto NDVI de la trama B5. Lo visualizo y lo guardo como un GeoTiff.

¿Cómo puedo crear algo así como un subconjunto de 800 x 600 píxeles de una imagen rasterizada TIFF? También quiero tomar ese TIFF y generar una imagen NDVI para ese subconjunto.

NOTA: Estoy trabajando con una imagen Landsat. La imagen tiene B5 al final del título del archivo.

Lo que he hecho hasta ahora:

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))

Esto funciona bien sin errores después de la ayuda de este sitio.

Quiero mostrar esto usando Spyder o Jupyter notebooks. Así que pensé en usar matplotlib e hice el código que fluye:

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

Haciendo esto se genera una ventana de matplotlib de 800x600, pero todo es de color púrpura, no estoy seguro de por qué se produce esto.

Ahora quiero poder mostrar esta imagen de 800x600. A continuación, quiero realizar un NDVI en ese subconjunto de imágenes de 800x600. A continuación, mostrar el subconjunto 800x600 imagen con NDVI mostrando.

Sé que la formuala lo es: NDVI = (NIR - rojo) / (NIR + rojo)

Pero, ¿cómo extraigo aquí el NIR y el rojo de esta única imagen Landsat?

Mi intento de eso:

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

Cuando ejecuto ese código para las bandas me sale el error:

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

Cuando ejecuto este código:

print(w.count)

Devuelve "1".

¿Esto significa que la imagen del Landsat sólo tiene una banda? Pero para hacer NDVI ¿no necesito 3 bandas?

Estoy pensando en escribir algún código como este para obtener el NDVI de ese raster. Pero no estoy seguro de cómo ir sobre la extracción de las bandas:

# 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)

Este código no funciona porque las bandas no están definidas como nada, así que no sé cómo definir las bandas para obtener el NDVI.

Después de esto no estoy seguro de cómo seguir con ambos mostrando y ahorrando la imagen renderizada.

Answer 1

Este es un ejemplo que utiliza lectura en ventana de rasterio para extraer un subconjunto de 800x600 de las bandas roja y NIR de una banda de ocho Imagen del Landsat 8, generando el NDVI y graficándolo en escala de grises.

from rasterio.windows import Window
import rasterio as rio

import  matplotlib.pyplot as plot
import numpy as np

with rio.open(r"LC80960722016224LGN00.tif") as src:
    red = src.read(4, window=Window(2500, 2500, 800, 600)).astype(np.float32)  # 2500, 2500 is roughly in the centre of my image
    nir = src.read(5, window=Window(2500, 2500, 800, 600)).astype(np.float32)
    red[red == src.nodata] = np.nan  # Convert NoData to NaN
    nir[nir == src.nodata] = np.nan  # Convert NoData to NaN

ndvi = (nir - red) / (nir + red)
vmin, vmax = np.nanpercentile(ndvi, (1,99))  # 1-99% contrast stretch

img_plt = plot.imshow(ndvi, cmap='gray', vmin=vmin, vmax=vmax)
plot.show()