Creating subset of TIFF image and generating its NDVI in Python?

Question

I have two Raster images, one from band4 with a B4 at the end and another with from band5 with B5 at the end. I want to subset the B5 raster to 800x600, then display it and save it as a GeoTiff. Then I want to compute the NDVI (I assume I'll need both the B4 and B5 to do this, but not sure). Then I want to display the NDVI subset of the B5 raster. Display it and save it as a GeoTiff.

How would I create something like a 800 x 600 pixel subset of a TIFF raster image? I want to also take that TIFF and generate an NDVI image for that subset.

NOTE: I am working with a Landsat image. The image has B5 at the end of the file title.

What I've done so far:

import rasterio
from rasterio.windows import Window 
import matplotlib.pyplot plt # for later use

with rasterio.open('MyRasterImage.tif') as src:
w = src.read(1, window=Window(0, 0, 800, 600))

This runs fine without error after help from this site.

I want to display this using Spyder or Jupyter notebooks. So I thought to use matplotlib and did the flowing code:

# Plot
plt.imshow(w)
plt.show()

Doing this generates a 800x600 matplotlib window, but it's all purple, not sure why its producing this.

Now I want to be able to display this 800x600 image. Then after that I want preform an NDVI on that subset 800x600 image. Then display the subset 800x600 image with NDVI showing.

I know the formuala is: NDVI = (NIR - red) / (NIR + red)

But how do I extract out NIR and red here from this single Landsat image?

My attempt at that:

band1 = dataset.read(1)
band2 = dataset.read(2)
band3 = dataset.read(3)
print(band[2])

When I run that code for the bands I get the error:

rasterio indexerror: band index 2 out of range (not in (1,))

When I run this code:

print(w.count)

It returns '1'.

So this means that the Landsat image only has one band? But in order to do NDVI don't I need 3 bands?

I am thinking of writing some code like this to get the NDVI from that raster. But not sure how to go about extracting out the bands:

# We handle the connections with "with"
with rasterio.open(bands[0]) as src:
    b3 = src.read(1)

with rasterio.open(bands[1]) as src:
    b4 = src.read(1)

# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

# Calculate NDVI
ndvi = (b4.astype(float) - b3.astype(float)) / (b4 + b3)

This code doesn't work because bands isn't define as anything so I don't know how to define bands to get the NDVI.

After this I am not sure how to go about both displaying and saving the rendered image.

Answer 1

Here's an example that uses rasterio windowed reading to extract an 800x600 subset of the red and NIR bands from an eight band Landsat 8 image, generating NDVI and ploting it in greyscale.

from rasterio.windows import Window
import rasterio as rio

import  matplotlib.pyplot as plot
import numpy as np

with rio.open(r"LC80960722016224LGN00.tif") as src:
    red = src.read(4, window=Window(2500, 2500, 800, 600)).astype(np.float32)  # 2500, 2500 is roughly in the centre of my image
    nir = src.read(5, window=Window(2500, 2500, 800, 600)).astype(np.float32)
    red[red == src.nodata] = np.nan  # Convert NoData to NaN
    nir[nir == src.nodata] = np.nan  # Convert NoData to NaN

ndvi = (nir - red) / (nir + red)
vmin, vmax = np.nanpercentile(ndvi, (1,99))  # 1-99% contrast stretch

img_plt = plot.imshow(ndvi, cmap='gray', vmin=vmin, vmax=vmax)
plot.show()