I am trying to write my first ArcPy script to iterate a Band Math process with requires a few steps. The uniterated, working script for a single group of bands in one Landsat scene is:

#import modules arcpy and string
import arcpy, string

#import the environment settings and spatial analyst extension
from arcpy import env
from arcpy.sa import *

#Check out the Spatial Analyst extension (must be available to check out!)

#Add your workspace here
env.workspace = r'C:\WORKSPACE'
env.overwriteOutput = True 

#Define inputs from TIFFs
Band2 = 'LANDSAT_B2.tif' #Band 2
Band4 = 'LANDSAT_B4.tif' #Band 4
Band5 = 'LANDSAT_B5.tif' #Band 5

#Band Ratios
R42 = arcpy.sa.Divide(Band4, Band2)
R52 = arcpy.sa.Divide(Band5, Band2)
R54 = arcpy.sa.Divide(Band5, Band4)
print 'Ratioed Bands'

a = 52
b = 11
c = 33.5
d = 120

#Algorithm processing
ALG = a - (b * R42) + (c * R52) - (d * R54)

#Saving final product
result = 'ALG.tif'
print 'Success'

My goal is to incorporate these steps into a for loop, so that I can run a script on a folder or set of folders rather than retyping each different scene's filename into the beginning of the script.

My plan is to separate the original bands, which are individual rasters, into folders and call them as raster objects:


and then begin by iterating the band ratio process through the folders, so that the first band 2 in the band 2 folder is ratioed with the first band 4 in the band 4 folder, and the second band 2 and 4 are ratioed, etc:

for idx in range(0, len(B2rasters) - 1):
    b2 = B2rasters[idx]
    b4 = B4rasters[idx]
output = b4 / b2
result = output + idx
print result

I don't feel confident in my ability to construct the iteration even after trawling through many examples of shorter algorithm iterations, as none of the examples I found use as many different rasters for multiple steps.

  • 2
    May I suggest this might be a bit too much to tackle for a first time user. Familiarize yourself with iteration first, like the example in resources.arcgis.com/en/help/main/10.2/index.html#//… and work up from there. I can see that your script just needs a for raster in arcpy.ListRasters loop but providing that without your understanding of why could lead to frustration when it doesn't work and you can't see why it isn't doing what its supposed to. – Michael Stimson Nov 21 '17 at 3:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.