ArcPy: Create NDVI with Raster Calculation and Map Algebra

Question

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?

Answer 1

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")))