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".

• is there a reason why you have to copy the bands to new rasters? It is more efficient to just use the original bands in your map algebra directly. Jun 8, 2015 at 9:56
• There's no reason for it. I'm just used to processing single band rasters. Jun 8, 2015 at 10:55
• If the "only" reason for band extraction is to apply your Python code above, there is a better way of doing this. You can iterate your rasters (assuming they are located in a/same folder) in Model Builder and calculate NDVI via Raster Calculator. Additionally your NDVI formula looks incorrect, please check the precedence of arithmetic operators. Jun 8, 2015 at 14:08

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
``````
• Finally, the above code has worked (after some tries), my next step is to batch implement this script. Jun 9, 2015 at 3:51

You can access an individual band from a multi-band raster by using the Make RasterLayer tool. This would create an "in-memory" raster which you can process with, it is not a new dataset.

• +1 Simpler than my answer if the band name is unknown. Jun 10, 2015 at 0:51