arcpy: For loop to sequence Processing of two items in list using ArcPy?

Question

I have multiple rasters, organized by years. I.e. raster y1 = from year 1 raster y2 = from year 2, etc...).

All of rasters contain values 0-11, representing forest (0-6) and disturbance (7-11). I would like to:

• Identify the forest on scene in y1 (by Reclassify tool)

• Use the forest from y1 to Extract by mask raster y2

• Identify the forest on scene in y2 (by Reclassify tool)

• Use the forest from y2 to Extract by mask raster y3

• Identify the forest on scene in y3 (by Reclassify tool)

• Use the forest from y3 to Extract by mask raster y4

etc for every year

However, I can't imagine how can I automate this process? I tried to run for loop and advance the raster index value by one to move to following raster in a list of rasters. Instead I have created unfinite loop, where the for loop always start with the first raster in a list... Or should I use another type of loop?

My approach:

# Import system modules
import arcpy
from arcpy.sa import *

# Set environment settings
wd = "C:/Users/Projects/2017_deforest"
env.workspace = wd + "/raw/SAO"
arcpy.env.overwriteOutput = True

# Set local variables
rasters = arcpy.ListRasters("sao*", "IMG")
reclassField = "Value"
remap = RemapRange([[0,0, "NODATA"],
                    [1,6,1],  # forest in that year
                    [7,11, "NODATA"], # disturbance 
                    [12,255, "NODATA"]])

# Check out the ArcGIS Spatial Analyst extension licence
arcpy.CheckOutExtension("Spatial")

# Execute Reclassify
for raster in rasters:

    # Identify forest in raster y1
    tempForest = Reclassify(raster, reclassField, remap, "NODATA")


     # Create list index to move to following raster in a list
    rasNumber = 1 
    inRaster2 = rasters[rasNumber]

    # Execute Extract by mask
    outName = wd + "/outputSAO/" + inRaster2
    outExtractByMask = ExtractByMask(inRaster2, tempForest)

    # advance the index by one to move to another element in  a list
    rasNumber = rasNumber + 1 

    # Save the output
     outExtractByMask.save(outName)

# Check in the ArcGIS Spatial Analyst extension license
arcpy.CheckInExtension("Spatial")

Answer 1

The problem is that you re-define

rasNumber = 1 

in each iteration, so you will be accessing the same raster in each iteration. This would not result in infinite loop though, because the for loop will exhaust once all raster files are iterated (of course you would get wrong results as the same - the first - raster was used). You should move this line outside of the for loop.

rasNumber = 1 

for raster in rasters:

When you will run your modified script, you would get an IndexError on the line inRaster2 = rasters[rasNumber] when you are processing the last file. This is because there is no raster that comes after it, so you cannot get from the list. This can be solved either by wrapping the indexing operation into a try/except statement:

try:
    inRaster2 = rasters[rasNumber]
except:
    break

Alternatively, you could just check the item being accessed:

if rasNumber > len(rasters):
    break

The for loop will break when you would start processing the last raster (you cannot extract by mask the non-existing raster).

Apart from that, some tips on the code that will make it safer:

• Don't use string concatenation as in outName = wd + "/outputSAO/" + inRaster2 and env.workspace = wd + "/raw/SAO". Instead, use the os.path.join to safely join folder paths and file names.

• Don't import all objects from the arcpy.sa module using from arcpy.sa import * because this takes extra time every time you run the script.

Answer 2

Thanks to @alextereshenkov answer I have put together my code. As it took me a while where exactly to put individual statements, here I attach the whole script for the record:

# Import system modules
import arcpy, os

# Set environment settings
wd = "C:/Users/Projects/2017_deforest"
arcpy.env.workspace = os.path.join(wd,"raw/SAO")
arcpy.env.overwriteOutput = True

# Set local variables
rasters = arcpy.ListRasters("sao*", "IMG")
reclassField = "Value"
remap = arcpy.sa.RemapRange([[0, 0, "NODATA"],  # need to specify the Spatial extension 
                    [1, 6, 2],  # forest in that year
                    [7, 11, "NODATA"], # disturbance 
                    [12, 255, "NODATA"]])


# Check out the ArcGIS Spatial Analyst extension licence
arcpy.CheckOutExtension("Spatial")

# define the index of the first raster in the list
rasNumber = 0

for raster in rasters:

    try:

        # Identify forest in raster y1
        tempForest = arcpy.sa.Reclassify(raster, reclassField, remap, "NODATA")

        # Execute Extract by mask        
       # print "(1) rasNumber is " + str(rasNumber)

        # Define the position of the second raster in the list
        inRaster2 = rasters[rasNumber]

        print "inRaster2 is " + str(inRaster2)

        outExtractByMask = ExtractByMask(inRaster2, tempForest)

        # Advance raster position
        rasNumber = rasNumber + 1
        print " (2)rasNumber is " + str(rasNumber)

        # Save the output
        outName = os.path.join(wd, "outputSAO", inRaster2)
        outExtractByMask.save(outName)

    except:

        # Break into code after overapassing the number of rasters in a list
        break  


# Check in the ArcGIS Spatial Analyst extension license
arcpy.CheckInExtension("Spatial")

print "DONE"