Calculates NDSSI from multispectral composite images from a folder

[MohsinMeraj]

I an new to arcpy and python. I have been trying to calculate NDSSI from 100 composite raster images from a folder using python script. Any how i tried my luck using a script but its not working. Taking list of rasters from one workspace and calculate NDSSI and store results in .tif format in another folder. 

import arcpy, string
from arcpy import envfrom arcpy.sa import*
arcpy.CheckOutExtension("spatial")
env.workspace = r'E:\Landsat All Images\Processing\LT5\Composite'
outws = r'E:\Landsat All Images\Processing\LT5\NDSSI'
rasters = arcpy.ListRasters("*.tif")
for raster in rasters:
Blue = raster + "\Band_1"NIR = raster + "\Band_4"Num = arcpy.sa.Float(Raster(Blue) - Raster(NIR))Denom = arcpy.sa.Float(Raster(Blue) + Raster(NIR))NIR_eq = arcpy.sa.Divide(Num, Denom) NIR_eq.save(outws)
print "Processing complete"

[Luke_Pinner]

Please never just say "it's not working". You are much more likely to get an answer by specifying what should happen and what actually happens, i.e does your script actually run at all, does it output the wrong result, does it output nothing, does it crash with an exception, what was the exception message, does python crash, does any smoke come out of your pc, does your pc vanish in a fiery explosion? Etc...

[DanPatterson_Retired]

You will need to format your code using the appropriate syntax highlighter.... It won't run as shown and impossible to dissect any errors that may be due to code structure or format.... have a look at

/blogs/dan_patterson/2016/08/14/script-formatting 

[MohsinMeraj]

Thank You Dan for Guidance....I am beginner in programming. Here is code with using online syntax highlighter.

# Calculates NDVI from multispectral imagery

import arcpy, string

from arcpy import env
from arcpy.sa import*

arcpy.CheckOutExtension("spatial")

env.workspace = r'E:\Landsat All Images\Processing\LT5\Composite'

outws = r'E:\Landsat All Images\Processing\LT5\NDSSI'

rasters = arcpy.ListRasters("*.tif")


for raster in rasters:

     # You may need to change the band combinations.  
     # This is for 4-band NAIP imagery or Landsat TM.
     Blue = raster + "\Band_1"
     NIR = raster + "\Band_4"

     Blue_out = "Blue.tif"
     NIR_out = "NIR.tif"

     arcpy.CopyRaster_management(Blue,Blue_out)
     arcpy.CopyRaster_management(NIR, NIR_out)

     Num = arcpy.sa.Float(Raster(Blue) - Raster(NIR))
     Denom = arcpy.sa.Float(Raster(Blue) + Raster(NIR))
     NIR_eq = arcpy.sa.Divide(Num, Denom)

     NIR_eq.save(outws)
print "Processing complete"‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

If memory serves, you will have to specify the filename and not just the environment workspace as in this helpfile example Raster—Help | ArcGIS Desktop 

Perhaps NIR_eq.save(outws+"/output.tif"

[XanderBakker]

Another thing to keep in mind, and what will probably create problems, is that fact that you save the Blue and NIR for each TIF to the same output name (in you input folder). Since you have not specified the "arcpy.env.overwriteOutput = True" it will produce an error for the second raster in the list.

[XanderBakker · Accepted solution]

I just executed the code below and it worked for me:

# Calculates NDVI from multispectral imagery
import arcpy
import os

from arcpy import env
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

# env.workspace = r'E:\Landsat All Images\Processing\LT5\Composite'
inws = r'C:\Esri EPM\Drone2Map\EPM\products\2D\ortho_tiles'
outws = r'C:\Esri EPM\Drone2Map\EPM\products\2D'
env.workspace = inws

rasters = arcpy.ListRasters("*.tif")
for raster in rasters:
    # You may need to change the band combinations.
    # This is for 4-band NAIP imagery or Landsat TM.
    Blue = arcpy.Raster(os.path.join(inws, raster + "\Band_1"))
    NIR = arcpy.Raster(os.path.join(inws, raster + "\Band_4"))

    NIR_eq = arcpy.sa.Float(Blue - NIR) / arcpy.sa.Float(Blue + NIR)
    out_name = raster.lower().replace('.tif', '_NDVI.tif')
    out_ras = os.path.join(outws, out_name)
    NIR_eq.save(out_ras)

print "Processing complete"‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍