Problema no cálculo de NDVI usando Rasterio Python

Question

Estou tentando calcular o NDVI usando duas imagens raster cortadas do Landsat 7 (NIR e bandas vermelhas cortadas usando o arquivo de máscara) usando o seguinte 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)

No código acima, ambas as bandas (Red e NIR) são de formas diferentes (linhas e colunas diferentes). Depois de executar o código acima, estou recebendo a mensagem "ValueError: operandos não puderam ser transmitidos junto com formas (2730,3635) (2731,3660)".

Mas quando o mesmo cálculo de NDVI que estou tentando fazer no ArcMap (usando a calculadora raster), o NDVI está sendo calculado.

Alguém pode me ajudar a resolver esse erro.

Answer 1

Para resolver seu problema, você precisa garantir que as grades cubram a mesma área e tenham as mesmas dimensões. Um método para conseguir isso é com o método reproject_match em rioxarray(extensão xarray geoespacial alimentada 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()