ArcPy: Create NDVI with Raster Calculation and Map Algebra

My target is to create NDVI Rasters but I'm stuck with python loops, since I'm new to that, so sorry for a maybe basic question.

I have 2 folders, one with rasters (tif) with the red band, the other folder with the same rasters with the infrared band.

Folder A: raster1.tif (red band) raster2.tif (red band) raster3.tif (red band)

Folder B: raster1.tif (infrared band) raster2.tif (infrared band) raster3.tif (infrared band)

I want to create NDVI rasters with the following formula: (NIR Band-Red Band)/(NIR Band+Red Band)

Of course, raster 1 from folder A has to match raster 1 from folder B.

I know I have to use map algebra (not Raster Calculator), but I'm already stuck with the loop, that the correct rasters are matching from folder A with folder B.

The following code doesn't produce the correct match. I tried it with a loop in a loop:

import arcpy
from arcpy.sa import*
arcpy.CheckOutExtension('Spatial')

arcpy.env.workspace = r"workspacetoredbandfolder"
liste_rotesband=arcpy.ListRasters()

arcpy.env.workspace = r"workspacetonirbandfolder"
liste_nirband=arcpy.ListRasters()

for rotesband in liste_rotesband:
print rotesband
for nirband in liste_nirband:
print nirband
break

but the result is:

raster 1- red band
raster 1- nir band
raster 2- red band
raster 1- nir band
raster 3- red band
raster 1- nir band

It should be:

raster 1- red band
raster 1- nir band
raster 2- red band
raster 2- nir band
raster 3- red band
raster 3- nir band

After this step, I was thinking about using the map Algebra to calculate my NDVI values.

Edit: I Tried it now with the following code:

import arcpy, os, math
arcpy.CheckOutExtension("Spatial")

workspaceRotesBand=r"workspacetoredbandfolder"
workspaceNIRBand=r"workspacetonirbandfolder"
OutputWS=r"workspacetooutput_ndvi_tiffs"

arcpy.env.workspace = workspaceRotesBand

rasterlist = arcpy.ListRasters()

for r in rasterlist:
r1 = arcpy.sa.Raster(r)
r2 = arcpy.sa.Raster(os.path.join(workspaceNIRBand,r))
result = Float((Raster("r2")-Raster("r1"))/(Raster("r2")+Raster("r1")))
result.save(os.path.join(OutputWS,r))

Now it seems that the linking works but I get the Error in the result line at the end: "'Float' is not defined". If I remove the Float I get the Error "'Raster' is not defined".

Is my map algebra expression correct?

If you want to get the matching rasters you can zip both lists so you iterate over a list of tuples containing two rasters (one red, one nir) at a time.

for rotesband, nirband in zip(liste_rotesband,liste_nirband):
print 'red band: {}'.format(rotesband)
print 'nir band: {}'.format(nirband)

Regarding your edit, you are attempting to use the Float() function. However, you are not importing it. You have two options:

1) Write from arcpy.sa import * on top of your script, or

2) Write arcpy.sa.Float((Raster("r2")-Raster("r1"))/(Raster("r2")+Raster("r1")))

• thank you @Marcelo Villa: after adding the line "from arcpy.sa import *", I am one step further but then the error: "000732: Input Raster: Dataset r2 does not exist or is not supported" appears. Any idea? – GIS7 Aug 13 at 14:04
• @GIS7 you already defined the rasters in r1 and r2. What happens if you do the following: Float((r2 - r1) / (r2+r1)) ? – Marcelo Villa Aug 13 at 14:21
• thank you Marcelo Villa! Now the tools doesnt produce the error. The result is still not what i want since i now get the values -1 and 0 as a result and not -1 till +1 with values in between (like 0.26). I just get 0 and -1. The result ist 32 bit floating point... – GIS7 Aug 14 at 6:09
• Now it works: Float((r2-r1))/Float((r2+r1)) Thank you very much @Marcelo Villa. – GIS7 Aug 14 at 6:41