I have two Raster images, one from band4 with a B4 at the end and another with from band5 with B5 at the end. I want to subset the B5 raster to 800x600, then display it and save it as a GeoTiff. Then I want to compute the NDVI (I assume I'll need both the B4 and B5 to do this, but not sure). Then I want to display the NDVI subset of the B5 raster. Display it and save it as a GeoTiff.
How would I create something like a 800 x 600 pixel subset of a TIFF raster image? I want to also take that TIFF and generate an NDVI image for that subset.
NOTE: I am working with a Landsat image. The image has B5 at the end of the file title.
What I've done so far:
import rasterio from rasterio.windows import Window import matplotlib.pyplot plt # for later use with rasterio.open('MyRasterImage.tif') as src: w = src.read(1, window=Window(0, 0, 800, 600))
This runs fine without error after help from this site.
I want to display this using Spyder or Jupyter notebooks. So I thought to use matplotlib and did the flowing code:
# Plot plt.imshow(w) plt.show()
Doing this generates a 800x600 matplotlib window, but it's all purple, not sure why its producing this.
Now I want to be able to display this 800x600 image. Then after that I want preform an NDVI on that subset 800x600 image. Then display the subset 800x600 image with NDVI showing.
I know the formuala is: NDVI = (NIR - red) / (NIR + red)
But how do I extract out NIR and red here from this single Landsat image?
My attempt at that:
band1 = dataset.read(1) band2 = dataset.read(2) band3 = dataset.read(3) print(band)
When I run that code for the bands I get the error:
rasterio indexerror: band index 2 out of range (not in (1,))
When I run this code:
It returns '1'.
So this means that the Landsat image only has one band? But in order to do NDVI don't I need 3 bands?
I am thinking of writing some code like this to get the NDVI from that raster. But not sure how to go about extracting out the bands:
# We handle the connections with "with" with rasterio.open(bands) as src: b3 = src.read(1) with rasterio.open(bands) as src: b4 = src.read(1) # Allow division by zero numpy.seterr(divide='ignore', invalid='ignore') # Calculate NDVI ndvi = (b4.astype(float) - b3.astype(float)) / (b4 + b3)
This code doesn't work because bands isn't define as anything so I don't know how to define bands to get the NDVI.
After this I am not sure how to go about both displaying and saving the rendered image.