pythonarraysnumpyrasterio

Calculating NDVI using python results in all zeros

I want to calculate the ndvi from a Sentinel-2 image.

import os
import numpy as np
import rasterio as rio

# suppress true divide warning numpy for 0 divide
np.seterr(divide='ignore', invalid='ignore')

red_f = absolute/path/to/band/4
nir_f = absolute/path/to/band/8

def calc_ndvi():
    
    with rio.open(red_f) as src:
        red = src.read()
        red = red.astype(np.float64)
    
    with rio.open(nir_f) as src:
        nir = src.read()
        nir = red.astype(np.float64)

    ndvi = np.divide((nir - red),(nir + red))
    
    return ndvi


ndvi = calc_ndvi()

The 'red' and 'nir' originally get loaded in as 'Array of uint16' with a shape of (1, 10980, 10980). I convert this to a float before the calculation using astype. As far as I know it's not necessary to flatten the array to a 2d shape. I have tried this but this didn't work.

The result unfortunately is an array completely filled with 0's.

What am I doing wrong?

Solution

  • You have a typo:

    nir = red.astype(np.float64)

    Should be:

    nir = nir.astype(np.float64)

    In:

    ndvi = np.divide((nir - red),(nir + red))

    You are really doing:

    ndvi = np.divide((red - red),(red + red))

    Which results in array of 0