I want to compute for an NDVI using rasterio. I must use two raster images that are adjacent to each other. They do not have the same sizes. They have a great overlap area which I can calculate NDVI on. One image is slightly shifted from the other. How to calculate raster that are not of the same sizes and does not overlap entirely?

Data: https://drive.google.com/drive/folders/1J3gmJVObv0LVoiUM_JhziD8ip19IB3l3?usp=sharing

EDIT 1: I have this to start with:

red_filename = "data/red.tif"
nir_filename = "data/nir.tif"

red_dataset = rasterio.open(red_filename)
nir_dataset = rasterio.open(nir_filename)

red_band = red_dataset.read(1)
nir_band = nir_dataset.read(1)

fig, (ax1, ax2) = pyplot.subplots(1,2, figsize=(10,12))
show((red_dataset, 1), ax=ax1, cmap='Reds', title='red_band')
show((nir_dataset, 1), ax=ax2, cmap='Blues', title='nir_band')
  • Are the two images acquired from the same sensor?
    – Aaron
    Nov 25 '20 at 16:58
  • Yes @Aaron. They are are adjacent captures of the same sensor.
    – Nikko
    Nov 25 '20 at 17:15
  • Could you please provide a screenshot highlighting the issue you are observing?
    – Aaron
    Nov 25 '20 at 17:23
  • I don't have an error. I need to know how to calculate the ndvi from the two raster. The two images does not have the same shapes. I want to be able to compute ndvi where the two images overlap.
    – Nikko
    Nov 25 '20 at 17:30
  • 1
    I think I understand your issue. Sounds like you have two images with separate bands (e.g. [/image1-band4.tif, /image1-band3.tif], [/image2-band4.tif, /image2-band3.tif]. If this is the case, you can indeed mosaic the two sets of images, so you will end up with, for example, a mosaiced band 4 and a mosaiced band 3. Then calculate NDVI. More details: automating-gis-processes.github.io/CSC18/lessons/L6/…
    – Aaron
    Nov 25 '20 at 17:51

I was able to do a stacking using osgeo/gdal. When stacked, it will generate a new image with a bounding box that extends to both image. After that you can read both band from the stacked image and perform the calculation. Each band from the stacked image will have 0 values on the extent which it not have original values.

from osgeo import gdal
red_filename = "data/red.tif"
nir_filename = "data/nir.tif"
outvrt = '/vsimem/stacked.vrt' #/vsimem is special in-memory virtual "directory"
outtif = 'data/stacked.tif'
tifs = [nir_filename, red_filename] 

outds = gdal.BuildVRT(outvrt, tifs, separate=True)
outds = gdal.Translate(outtif, outds)

stack_filename = "data/stacked.tif"
stack_dataset = rasterio.open(stack_filename)

nir_band = stack_dataset.read(1)
red_band = stack_dataset.read(2)

fig, (ax1, ax2) = pyplot.subplots(1,2, figsize=(10,12))
show((red_band), ax=ax1, cmap='Reds', title='red_band')
show((nir_band), ax=ax2, cmap='Blues', title='nir_band')

ndvi = (nir_band - red_band) / (nir_band + red_band)
fig, ax = pyplot.subplots(figsize=(6,8))
show((ndvi), ax=ax, cmap='RdYlGn', title='NDVI')

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.