0

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)
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)
  • Thank you. The files I am working with I have downloaded on to my desktop. I am going about it different than the tutorial. I'd just like to compute NDVI with the two TIFFS I have on my desktop. Since those don't have overviews. Do you know of another way of still calculating the NDVI without them? – yuen2 Nov 19 '18 at 4:23
  • Sure, I just updated the answer to include how to calculate NDVI from the files you already have downloaded. – James Nov 19 '18 at 22:11
  • Thank you! Much appreciated! I am able to display the NDVI fine I named it naip_ndvi. How would I save it? I made an attempt at the end of my code. Every time I try to save it gives me a black screen. – yuen2 Nov 19 '18 at 22:18
  • Based on the source code itself, I think you've saved the file correctly! Most of the time, though, you'll want to visualize it in a GIS application such as QGIS (qgis.org). Viewing your raster in your OS's image viewer application will usually look like a black rectangle. – James Nov 20 '18 at 5:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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