python - Access individual bands and use them in map algebra

I have a list of multiband rasters, consisting of: Red=Band 1, NIR=Band 2, SWIR1=Band 3 and SWIR2= Band 4. What I want is to access only the red and NIR bands and save it to new file. Then, I'll be using them my Map Algebra calculations below.

(Note that accessing and saving the Red and NIR bands portion is still missing from my code)

    import arcpy
    from arcpy import env

    from arcpy.sa import *
    env.overwriteOutput = True

    #Set the current workspace, assuming that all red and nir bands are already 
    #extracted and saved as individual bands
    env.workspace = (r"C:\thesis\hansenwipfiles\first and last v1.0")

    #Raster List (Assuming the rasters here are already single bands)
    RedbandList = ["red20N_120E.tif", "red10N_120E.tif", "red20N_110E.tif",     
                   "red10N_110E.tif"]  


    NIRbandList = ["NIR20N_120E.tif", "NIR10N_120E.tif", "NIR20N_110E.tif",     
                   "NIR10N_110E.tif"]

    for range in (0,4):
        ndvi = Raster(NIRbandList[f])-Raster(RedbandList[f])/\
                  Raster(NIRbandList[f])+Raster(RedbandList[f])
        ndviFloat = float(ndvi)
        ouputName = ndviFloat
        ndviFloat.save(ouputName)


    print "Finish!"

I have not tested the code if it works, though no errors found after running the "Test Module".

Answer

You can access individual bands by joining the band name to the raster path - i.e. path/to/raster/band_name.

Often using path/to/raster/Band_[band number] works (i.e os.path.join(rasterpath, 'Band_1'), but not always. I use Landsat 8 imagery quite a bit and ArcGIS names the bands 'CoastalAerosol', 'Blue', 'Green', 'Red', etc...

If you don't know what the band names are, you can check the raster Properties->KeyMetadata tab -> Source Band Index (note use ArcCatalog or the Catalog window in ArcMap to check properties, don't do it from the ArcMap TOC).

Alternatively, you can get the band names in code by setting the current workspace to the multiband raster, then calling arcpy.ListRasters() to get a list of the band names.

For example:

import os
import arcview, arcpy
from arcpy.sa import Raster, Float

arcpy.CheckOutExtension('spatial')

ls8ms = r'P:\temp\test.tif' #Some landsat 8 multispectral data

#Set workspace to multiband raster, then list rasters to get band names
#There must be a more obvious way of doing this...

arcpy.env.workspace = ls8ms
bands = [Raster(os.path.join(ls8ms, b)) for b in arcpy.ListRasters()]

# Set workspace somewhere sensible
arcpy.env.workspace = r'P:\temp'

# Do the calculation. 
# Note band numbers are 0 indexed. 
# Landsat 8 - Band 5 (NIR) = bands[4],  Band 4 (Red) = bands[3]
# Landsat 7 - Band 4 (NIR) = bands[3],  Band 3 (Red) = bands[2]

# For your data -  Band 2 (NIR) = bands[1],  Band 1 (Red) = bands[0]
ndvi = (Float(bands[4]) - bands[3]) / (Float(bands[4]) + bands[3]) #Landsat 8
# ndvi = (Float(bands[1]) - bands[0]) / (Float(bands[1]) + bands[0]) #Your data
ndvi.save(r'P:\temp\testndvi.tif')

arcpy.CheckInExtension('spatial')

You could wrap that up into a little function:

def get_bands(path_to_raster):
    """ Get a list of bands as Raster objects from a multiband raster """


    oldws = arcpy.env.workspace #Save previous workspace

    #Get raster objects from band names
    arcpy.env.workspace = path_to_raster
    bands = [Raster(os.path.join(path_to_raster, b)) for b in arcpy.ListRasters()]

    #Restore previous workspace
    arcpy.env.workspace = oldws


    return bands