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(2000,2000,800,600)

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

plt.figure(figsize=(6,8.5))
plt.imshow(subset)
plt.title(f'Band 5 Subset')
plt.xlabel('Column #')
plt.ylabel('Row #')





rast = rasterio.open(nirband)
rast2 = rasterio.open(redband)
nir = rast.read(1)
red = rast2.read(1)


red = red.astype(float)
nir = nir.astype(float)
np.seterr(divide='ignore', invalid='ignore')


ndvi = np.empty(rast.shape, dtype=rasterio.float32)
check = np.logical_or ( red > 0, nir > 0 )
ndvi = np.where ( check,  (1.0*(nir - red )) / (1.0*( nir + red )),-2 )


plt.figure(figsize=(6,8.5))
plt.imshow(ndvi)
plt.title(f'NDVI')
plt.xlabel('Column #')
plt.ylabel('Row #')