Using Rasterio, Matplotlib, and Numpy trying to calculate and plot NDVI of a raster image getting error: IndexError: list index out of range?

Question

First I create two variables and store my raster images in them and I use MatPlotLib to plot a subset of the original "nirband" raster image. That works fine. I am getting an IndexError: list index out of range error now. I am basing my code off the following tutorial found on GitHub:

https://github.com/geohackweek/tutorial_contents/blob/ba5e9443137a9aca87cdcdcd70e9e6a237cc64ba/raster/notebooks/rasterio-landsat-aws.ipynb

Here is the first section of the code which plots fine without error:

import rasterio
import matplotlib.pyplot as plt
import numpy as np

nirband = r"LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF"

redband =r"LC08_L1TP_015033_20170822_20170912_01_T1_B4.TIF"

#rasterio.windows.Window(col_off, row_off, width, height)
window = rasterio.windows.Window(1024, 1024, 800, 600)

with rasterio.open(nirband) as src:
    subset = src.read(1, window=window)

plt.figure(figsize=(15,8))
plt.imshow(subset)
plt.title(f'Band 5 Subset\n{window}')
plt.xlabel('Column #')
plt.ylabel('Row #')

All the code before this point runs fine, I was able to plot the above code using MatPlotLib just fine. It's the following code which produces an error message where I am trying to compute and display the NDVI of that raster subset image:

with rasterio.open(redband) as src:
    profile = src.profile
    oviews = src.overviews(1) # list of overviews from biggest to smallest
    oview = oviews[1]  # Use second-highest resolution overview
    print('Decimation factor= {}'.format(oview))
    red = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))

plt.imshow(red)
plt.colorbar()
plt.title('{}\nRed {}'.format(redband, red.shape))
plt.xlabel('Column #')
plt.ylabel('Row #')




with rasterio.open(nirband) as src:
    oviews = src.overviews(1) # list of overviews from biggest to smallest
    oview = oviews[1]  # Use second-highest resolution overview
    nir = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))

plt.imshow(nir)
plt.colorbar()
plt.title('{}\nNIR {}'.format(nirband, nir.shape))
plt.xlabel('Column #')
plt.ylabel('Row #')





def calc_ndvi(nir,red):
    '''Calculate NDVI from integer arrays'''
    nir = nir.astype('f4')
    red = red.astype('f4')
    ndvi = (nir - red) / (nir + red)
    return ndvi




ndvi = calc_ndvi(nir,red)
plt.imshow(ndvi, cmap='RdYlGn')
plt.colorbar()
# maybe edit it
plt.title('NDVI')
plt.xlabel('Column #')
plt.ylabel('Row #')

Here is the complete traceback with error message:

runfile('C:/Users/new/Desktop/RasterNDVI.py', wdir='C:/Users/new/Desktop')
Traceback (most recent call last):

  File "<ipython-input-15-5d04fa0ce75f>", line 1, in <module>
    runfile('C:/Users/new/Desktop/RasterNDVI.py', wdir='C:/Users/new/Desktop')

  File "C:\Users\new\Miniconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 668, in runfile
    execfile(filename, namespace)

  File "C:\Users\new\Miniconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/new/Desktop/RasterNDVI.py", line 27, in <module>
    oview = oviews[1]  # Use second-highest resolution overview

IndexError: list index out of range

Attempt at saving NDVI as a GeoTIFF file:

with rasterio.open("LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF") as src:
    naip_data_ras = src.read()
    naip_meta = src.profile



with rasterio.open('MyExample.tif', 'w',**naip_meta) as dst:
    dst.write(naip_ndvi, window=window)

Answer 1

Is it possible that you downloaded the source GeoTiff (LC08_L1TP_015033_20170822_20170912_01_T1_B4.TIF) but not the overviews? If you take a look at the index for this landsat scene, each landsat band has an associated .ovr file, which contains the overviews. If you download the .ovr file for the NIR band (band 5) and try re-running the program, you should be able to load the overviews as expected.

By the by, one of the things this tutorial is trying to demonstrate is the features within the Cloud-Optimized GeoTiff, which means that the example allows you to open a raster by its location on the web rather than having to download a file before opening it with rasterio. So, a complete COG-style opening of the raster would look like:

url = 'https://landsat-pds.s3.amazonaws.com/c1/L8/042/034/LC08_L1TP_042034_20170616_20170629_01_T1/LC08_L1TP_042034_20170616_20170629_01_T1_B4.TIF'
with rasterio.open(url) as src:
    oviews = src.overviews(1) # list of overviews from biggest to smallest
    oview = oviews[1]  # Use second-highest resolution overview
    nir = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))

Edit:

To calculate NDVI on the full raster that you already have downloaded, just read in the pixel values directly for each band:

nirband = r"LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF"
redband = r"LC08_L1TP_015033_20170822_20170912_01_T1_B4.TIF"
red = rasterio.open(redband).read(1)
nir = rasterio.open(nirband).read(1)

calc_ndvi(nir, red)