Problema nel calcolo dell'NDVI utilizzando Rasterio Python

Question

Sto cercando di calcolare l'NDVI usando due immagini raster ritagliate di Landsat 7 (NIR e bande rosse ritagliate usando un file maschera) usando il seguente codice:

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)

Nel codice sopra entrambe le bande (Rosso e NIR) sono di forme diverse (righe e colonne diverse). Dopo aver eseguito il codice sopra, ricevo il messaggio "ValueError: gli operandi non possono essere trasmessi insieme alle forme (2730,3635) (2731,3660)".

Ma quando lo stesso calcolo NDVI che sto cercando di fare in ArcMap (usando Raster Calculator), viene calcolato l'NDVI.

Qualcuno può aiutarmi a risolvere questo errore.

Answer 1

Per risolvere il problema, è necessario assicurarsi che le griglie coprano la stessa area e abbiano le stesse dimensioni. Un metodo per ottenere ciò è con il metodo reproject_match in rioxarray(estensione geospaziale xarray alimentata da 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()