Calculates NDSSI from multispectral composite images from a folder

1940
6
Jump to solution
12-31-2016 11:35 AM
MohsinMeraj
New Contributor

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"

 

0 Kudos
1 Solution

Accepted Solutions
XanderBakker
Esri Esteemed Contributor

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"

View solution in original post

6 Replies
Luke_Pinner
MVP Regular Contributor

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
MVP Emeritus

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
New Contributor

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"
0 Kudos
DanPatterson_Retired
MVP Emeritus

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
Esri Esteemed Contributor

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
Esri Esteemed Contributor

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"