Calculate NDVI using ArcPy Spatial Analyst

Question

My goal is to calculate NDVI using map algebra syntax.

import arcpy
from arcpy.sa import *

red_band_raster = (r"C:\B3.tif")
nir_band_raster = (r"C:\B4.tif")

#Calculating NDVI float raster
ndvi_raster = Divide(
    Float(Raster(nir_band_raster) - Raster(red_band_raster)), Float(Raster(nir_band_raster) + Raster(red_band_raster)))

#Saving raster:
ndvi_raster.save(r"C:\NDVI.tif")

Unfortunately, instead of varying NDVI values I get constant raster filled with "1.000000" and NoData

Answer 1

Are the red_band_raster and the nir_band_raster values integers? It seems likely that you should move the float() commands around each raster variable call in your equations, and not around the result of the addition/subtraction as you have it now.

Answer 2

I presume you're using ArcGIS Desktop which uses python 2x, not ArcGIS Pro which uses python 3x. Side note: this is one example of why it's important to specify the software you are using.

In python 2x, when both operands are integers, division returns an integer result rounded towards minus infinity. This was changed in Python 3x so division returns a float unless you force integer division using the // operator.

Here's a worked example to demonstrate:

import arcpy
from arcpy.sa import Raster, Int, Float, CreateConstantRaster
arcpy.CheckOutExtension('spatial')

red = Int(CreateConstantRaster(500))
nir = Int(CreateConstantRaster(1000))

#Or load from file
#Note load files into Raster objects just once,
#instead of each time you use them as you were
#in your original code

#red = Raster(r'path\to\red.tif')
#nir = Raster(r'path\to\nir.tif')

#Calculate NDVI using the python division operator instead
#of the ArcGIS "Divide" function so code is easier to read

#Integer division
ndvi = (nir - red) / (nir + red)

print ndvi.pixelType, ndvi.maximum, ndvi.minimum
#prints S32 (signed 32bit int) 0.0 0.0

#One operand floating point
ndvi = Float((nir - red)) / (nir + red)

print ndvi.pixelType, ndvi.maximum, ndvi.minimum
#prints F32 (32bit float) 0.333333343267 0.333333343267

Alternatively, you can force floating point division by using from __future__ import division at the very top of your script:

from __future__ import division

#imports etc...

red = Raster(r'path\to\red.tif')
nir = Raster(r'path\to\nir.tif')

#Integer division
ndvi = (nir - red) / (nir + red)

print ndvi.pixelType, ndvi.maximum, ndvi.minimum
#prints F32 (32bit float) 0.333333343267 0.333333343267