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

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.