Using rasterio and matplotlib how can I create a window for an NDVI image to fit into?

Question

I am using Rasterio, Numpy, and Matplotlib. I created a window with the dimensions of 800x600 and plotted the band 5 for a Landsat 8 image, that worked just fine. I also have the nir image so then I calculated the ndvi of the red and nir bands from the Landsat8 image. What I am trying to do is put my ndvi image inside a 800x600 window. How might this be done? So far or the ndvi image it's giving me default dimensions of 7000x7000, I want to create the or use the same window of Window(2000,2000,800,6000)

Here is my code so far:

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 #')

This creates two plots one which is 800x600, the other ranges from 7000x7000, I think this is the default for the image. I am trying to have the same dimensions for my ndvi image as I have for my nir band image.

Here is my attempt at saving the 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_band(1, naip_ndvi, window=window)

Answer 1

Try this code. Just read the first two bands in, windowed as you want them, then perform the ndvi calculation on those rasters, resulting in an ndvi image of the same window size.

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

with rasterio.open(redband) as src:
    red= src.read(1, window=window)

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

ndvi = np.empty(nir.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 )

Your problem was that you were re-loading your full red and nir bands as (7000 x 7000) rasters, then performing your ndvi calculation. What you should do (and what I've shown above) is to load in your red and nir bands for only the windowed area, then perform the ndvi calculation.