Problema al calcular NDVI usando Rasterio Python

Question

Estoy tratando de calcular NDVI usando dos imágenes rasterizadas recortadas de Landsat 7 (NIR y bandas rojas recortadas usando un archivo de máscara) usando el siguiente código:

import rasterio as rio
import numpy as np
import matplotlib.pyplot as plt

with rio.open(r'D:\clip_test_b3.tif') as src:
    red = src.read(1)   # (Rows, Columns) = (2731, 3660)
    
with rio.open(r'D:\clip_test_b4.tif') as src:
    nir = src.read(1)   # (Rows, Columns) = (2730, 3635)
    
np.seterr(divide = 'ignore', invalid = 'ignore')

ndvi = (nir.astype(float) - red.astype(float))/(nir + red)

plt.imshow(ndvi)

En el código anterior, ambas bandas (rojo y NIR) tienen diferentes formas (diferentes filas y columnas). Después de ejecutar el código anterior, recibo el mensaje "ValueError: los operandos no se pudieron transmitir junto con las formas (2730,3635) (2731,3660)".

Pero cuando estoy tratando de hacer el mismo cálculo de NDVI en ArcMap (usando la Calculadora de ráster), entonces se está calculando NDVI.

¿Alguien puede ayudarme a resolver este error?

Answer 1

Para resolver su problema, debe asegurarse de que las rejillas cubran la misma área y tengan las mismas dimensiones. Un método para lograr esto es con el método reproject_match en rioxarray(extensión geoespacial xarray impulsada por rasterio).

import rioxarray

red = rioxarray.open_rasterio("D:\clip_test_b3.tif")
nir_original = rioxarray.open_rasterio("D:\clip_test_b4.tif")

nir = nir_original.rio.reproject_match(red)
ndvi = (nir.astype(float) - red.astype(float))/(nir + red)

ndvi.plot()