1

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

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.

  • Thank you! How would I do it the first way? i'd like to have either another window or the same window and just have it set up like Window(2000,2000,800,600) and have the ndvi tucked in there. – yuen2 Nov 19 '18 at 19:18
  • Also, all I have access to right now is the two bands B5 and B4. So just using those I'm trying to create the same dimension box for NIR and NDVI. – yuen2 Nov 19 '18 at 19:19
  • Here is an image of the output from Spyder, I'd like to have both windows with the same dimensions and subset area, so col_off and col_on_ would have to be the same: imgur.com/a/gf0ohAp – yuen2 Nov 19 '18 at 19:29
  • @yuen2 Did you compute the NDVI yourself or get it from another location? Does it have the same dimensions as the Landsat images? – Jon Nov 19 '18 at 19:31
  • Yes NDVI is calculate from this line in the code down: ndvi = np.empty(rast.shape, dtype=rasterio.float32)…. now just trying to put that into a dimension box with Window(2000,2000,800,600) – yuen2 Nov 19 '18 at 19:35

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.