Issue in calculating NDVI using Rasterio Python

Question

I am trying to calculate NDVI using two clipped raster images of Landsat 7 (NIR & Red Bands clipped using mask file) using the following code:

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)

In the above code both the bands (Red & NIR) are of different shapes (different rows and columns). After running the above code I am getting message "ValueError: operands could not be broadcast together with shapes (2730,3635) (2731,3660) ".

But when same NDVI calculation I am trying to do in ArcMap (using Raster Calculator), then NDVI is getting calculated.

Can someone please help me out in solving out this error.

Answer 1

In order to solve your problem, you need to ensure the grids cover the same area and have the same dimensions. One method to achieve this is with the reproject_match method in rioxarray (geospatial xarray extension powered by 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()