Calculating NDVI using Rasterio, Python

Question

Below is my code for calculating NDVI using Landsat 8. I have used band 4 and 5, but I am not getting any output. How should I get the output image?

import numpy as np 
import rasterio
import subprocess
b4 = rasterio.open('C:\\Users\\ARSH\\Downloads\\LC08_L1TP_143051_20191219_20191226_01_T1\\LC08_L1TP_143051_20191219_20191226_01_T1_B4.tif')
b5 = rasterio.open('C:\\Users\\ARSH\\Downloads\\LC08_L1TP_143051_20191219_20191226_01_T1\\LC08_L1TP_143051_20191219_20191226_01_T1_B5.tif')
red = b4.read()
nir = b5.read()
np.seterr(divide='ignore', invalid='ignore')
# Calculate ndvi
ndvi = (nir.astype(float)-red.astype(float))/(nir+red)

# Write the NDVI image
meta = b4.meta
meta.update(driver='GTiff')
meta.update(dtype=rasterio.float32)

with rasterio.open('NDVI.tif', 'w', **meta) as dst:
    dst.write(ndvi.astype(rasterio.float32)

Answer 1

I adapted your code for my own raster layers (Landsat 5):

import numpy as np 
import rasterio

b4 = rasterio.open('/home/zeito/pyqgis_data/b3_ref.tif')
b5 = rasterio.open('/home/zeito/pyqgis_data/b4_ref.tif')
red = b4.read()
nir = b5.read()
np.seterr(divide='ignore', invalid='ignore')
# Calculate ndvi
ndvi = (nir.astype(float)-red.astype(float))/(nir+red)

# Write the NDVI image
meta = b4.meta
meta.update(driver='GTiff')
meta.update(dtype=rasterio.float32)

with rasterio.open('/home/zeito/pyqgis_data/NDVI.tif', 'w', **meta) as dst:
    dst.write(ndvi.astype(rasterio.float32))

and, after running it, I found out that the issue was due a missing parentheses in last instruction. Result can be visualized in following image:

NDVI values vary (-0.614977, 0,860435) as expected.

In your case (as Windows user), you should get the output image in:

C:\\Users\\ARSH\\

If you have some issues for founding it, put the full pass in following instruction:

.
.
.
with rasterio.open('NDVI.tif', 'w', **meta) as dst:
.
.
.