3

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. – user2856 Jun 8 '15 at 9:56
  • There's no reason for it. I'm just used to processing single band rasters. – brentiemapper Jun 8 '15 at 10:55
  • 1
    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. – fatih_dur Jun 8 '15 at 14:08
9

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. – brentiemapper Jun 9 '15 at 3:51
2

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. – user2856 Jun 10 '15 at 0:51

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.