1

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?

3

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
  • 1
    @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

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.