Calling individual bands from multiband.tif in python perform action and save output in mutilband.tif?

Question

I am trying to write a script that would go through my directory and open each *.tif file, evaluate each of their 120 bands and write the results out to another 120 multiband *.tif file. The only way I could figure out how to do this was to iterate through the input directory,

open each raster,

iterate through its bands,

evaluate each band for the threshold value (2),

if condition is met assign a number if condition is not met assign a value of zero,

save each band to an output.tif file

then use the composite function to collect all 120 output.tifs

But I'm stuck. When I try to iterate through the bands I get an error message (see below).

First is there an easier way to do this without creating 120 files for each of my 350 multiband input.tifs? If not, can you help me figure out how to iterate through each band?

This is the error message:

========= RESTART: C:/Python27/ArcGIS10.5/Scripts/StartJD.py ========= about to start cycling through rasters in the raster cycle at count 0

Traceback (most recent call last): File "C:/Python27/ArcGIS10.5/Scripts/sapflow_StartJD.py", line 36, in sapRas = Raster(inRas + "\Band_" + count) #get individual bands for input raster File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4329, in Plus in_raster_or_constant2) File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Utils.py", line 53, in swapper result = wrapper(*args, **kwargs) File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4326, in Wrapper ["Plus", in_raster_or_constant1, in_raster_or_constant2]) RuntimeError: ERROR 000732: Input Raster: Dataset \Band_ does not exist or is not supported

This is an example of the input file names within the directory, 26_sapsum_2011_converted (note no extension is listed but the accompanying file does have the extension 26_sapsum_2011_converted.tif.aux).

Here's my script:

# objective of this code is to assign the Julian day when sapsum = 2 and 
write a multiband output file for each year 
# Import arcpy module
import arcpy, os, sys
#import the environment settings and spatial analyst extension
from arcpy import env
arcpy.env.extent = "MAXOF"
# check out spatial analyst extension
if arcpy.CheckExtension('Spatial') == 'Available':
    arcpy.CheckOutExtension('Spatial')
else:
    arcpy.AddMessage("Error: Couldn't get Spatial Analyst extenstion, 
exiting")
sys.exit(1) 

from arcpy.sa import *

# Set the input/output directories
env.workspace = r"C:\Netcfd\BCSD\sapsum\sapsumtiff\Year\26_Y_20s"

# set number of bands
number_of_bands = 120
# Get a list of Rasters
rasters = arcpy.ListRasters('*.tif')

print 'about to start cycling through rasters'

for raster in rasters:
  # cycle through rasters in directory     
 inRas = Raster(raster) #open and rename raster input
 count = 0
 print "in the raster cycle at count", count
 for band_number in range(0, number_of_bands):
      count = count + 1
      sapRas = Raster(inRas + "\Band_" + count) #get individual bands for 
      # input raster
      print 'in the band number cycle and the current band number is', count
      output = con((sapRas == 2),count,0) #if sapRas is equal to 2 assign 
      #count if SapRas not equal to 2 assign a value of 0
      #save output in temporary file
      Rasname, Rasext = os.path.splitext(inRas)# sparate the name from the 
        #extension
      save.output(env.workspace + "/JD" + Rasname + "_comp.tif")
      for i in range(0, number_of_bands):
           listcomp = [] #create a list of output tifs
           for ras in arcpy.ListRasters('_comp.tif'):
                listcomp.append(ras)
           arcpy.CompositeBands_management(listcomp, env.workspace + 
          "Composite" + inRas + '.tif') #collect output tifs into a 
                                        #composite
print "Done! and look it's beer o'clock"

raw_input("Press <enter>") # this keeps the window open until you press 
                            "enter"

Answer 1

There's a few burning problems with your code.. Bands are 1 based, and it's better to keep one iterator rather than trying to use two, I suggest loop for band_number in range(1, number_of_bands + 1): and use band_number as the iterator.

As much as you're sure there's 120 bands if you use the Describe object properties to get the bandCount for looping you can use this script with more or less bands without changing... a good habit to get into is to make your code as generic as possible without overly complicating.

The syntax for saving a raster object is raster.save('full path') so save.output() isn't valid.

Start with an empty list for each input raster and append each raster in turn then composite bands at the end of the band_number loop.

I've changed the print statements to arcpy.AddMessage() just in case you want to use this from the toolbox - if you get into the habit now you won't have to change all the print statements later

Here's the guts of your loop:

for raster in rasters:
    # cycle through rasters in directory     
    d = arcpy.Describe(raster)
    number_of_bands = d.bandCount # the number of bands from the describe object
    arcpy.AddMessage("Processing raster {} for {} bands".format(raster,number_of_bands))
    listcomp = [] #create a list of output tifs, empty for each input raster

    for band_number in range(1, number_of_bands + 1):
            sapRas = arcpy.Raster("{}\\Band_{}".format(raster,band_number)) #get individual bands for input raster
            arcpy.AddMessage('in the band number cycle and the current band number is {}'.format( band_number))
            output = arcpy.sa.Con((sapRas == 2),band_number,0)  #if sapRas is equal to 2 assign count if SapRas 
                                                                                                                    #not equal to 2 assign a value of 0
            #save output in temporary file and add it to the list
            Rasname, Rasext = os.path.splitext(raster)# sparate the name from the extension
            output.save(os.path.join(env.workspace, "JD{}comp.tif".format(Rasname)))
            listcomp.append(os.path.join(env.workspace, "JD{}comp.tif".format(Rasname)))

    #collect output tifs into a composite
    arcpy.CompositeBands_management(listcomp, os.path.join(env.workspace,"Composite{}.tif".format(raster)))