2

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

  • 1
    @traggatmot, you were right, wrapping each raster in Float() solved the problem. Could you post your answer as an answer, or should I just post my fixed code? – Basile Aug 3 '17 at 15:01
4

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.

  • Thanks, turning inputs into Float(Raster(%input_tff%)) solved the problem – Basile Aug 4 '17 at 8:57
  • You may want to take some of Luke's advice (re: loading the rasters in once) and also mark his answer as 'the' answer - its more detailed. – traggatmot Aug 7 '17 at 13:22
5

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

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.