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.

• user2856: Can you edit your question to include your clip code, you're getting different sized outputs. And also rio info or gdalinfo outputs for the pre-clip dataset/s.

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

• user2856: +1 this answers the question and addresses the symptoms (and I wasn't aware of this useful functionality, so thanks), but I think if something has gone awry in the OPs clipping that they probably need to sort it out instead of just dealing with the different shapes afterwards.
• snowman2: True. You can use rioxarray to clip if you want to start from the beginning: gis.stackexchange.com/a/354798/144357 or corteva.github.io/rioxarray/stable/examples/clip_geom.html
• user2856: Thanks, I must have a look at rioxarray sometime.
• RRSC NGP: @snowman2. Thank for your help. It has solved my issue. I think there wasn't any issue with clipping. As the bands (Red & NIR) corresponds to the same "Path & Row" of Landsat-7, after making the polygon of the valid data for both the bands I found that both the polygons doesn't overlap exactly over each other. Because of this when I am masking the bands by using my study area the shapes of both the bands are not exactly same. Small edit in your solution "nir = nir_original.reproject_match(red)" should be changed to "nir = nir_original.rio.reproject_match(red)"
• snowman2: Glad to hear that it worked - and good catch with the missing rio. I just updated the answer.