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

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)

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:


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.

  • @Luke Thank you! I used the Reader code snippet, but am getting this error: w = src.read(1, window = Window(0, 0, 800, 600)) NameError: name 'Window' is not defined
    – yre
    Commented Nov 16, 2018 at 17:08
  • Thank you. I did: import rasterio THEN from rasterio import windows, from there it gives me the error: NameError: name 'Window' is not defined
    – yre
    Commented Nov 16, 2018 at 21:09
  • @Luke Thanks you, yes that worked! I did w = src.read(1, window = Window( 0,0,800,600 )). How do I display it and then save it for a file?
    – yre
    Commented Nov 17, 2018 at 14:36
  • @Luke I added details to my original question.
    – yre
    Commented Nov 17, 2018 at 16:24
  • gis.stackexchange.com/help/someone-answers
    – user2856
    Commented May 16, 2019 at 5:47

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

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.