pro calculate_savi
; Set compile options
compile_opt IDL2
; General error handler
CATCH, err
if (err ne 0) then begin
CATCH, /CANCEL
if OBJ_VALID(e) then $
e.ReportError, 'Please provide a valid multispectral file', /INFORMATION
MESSAGE, /RESET
return
endif
; Get ENVI session
e = ENVI(/CURRENT)
; Get the User Interface (UI) object
UI = e.UI
; Open ENVI's Quickbird Dataset
file = FILEPATH('qb_boulder_msi', ROOT_DIR=e.ROOT_DIR, SUBDIRECTORY = ['data'])
raster = e.OpenRaster(file)
; Now display the chosen raster
view= e.GetView()
layer = view.CreateLayer(raster, bands=[3,2,1])
; Return a file ID to pass to ENVI DOIT routines
fid = ENVIRasterToFID(raster)
; use FID to query file to get parameters needed for DOIT routines
ENVI_FILE_QUERY, fid, ns=ns, nl=nl, nb=nb, dims=dims, sname=sname
pos=lindgen(nb)
; set output directory
outDir = 'c:\testing\dataout\'
; calibrate to reflectance using Quac
ENVI_DOIT, 'ENVI_QUAC_DOIT', fid=fid, dims=dims, pos=pos, quac_sensor='QuickBird', $
out_name=outDir+sname+'_refl.dat', r_fid=refl_fid
; set up parameters for MATH_DOIT to calculate SAVI
; SAVI expression
exp = '((float(b4) - float(b3)) / ((float(b4) + float(b3)) + 0.5)) * 1.5'
in_pos = [2,3]
in_fid = [refl_fid, refl_fid]
; call MATH_DOIT
ENVI_DOIT, 'MATH_DOIT', fid=in_fid, pos=in_pos, dims=dims, exp=exp, $
out_name=outDir+sname+'_SAVI.dat', r_fid=savi_fid
;get raster object from FID
newRaster = ENVIFidToRaster(savi_fid)
; get all the layers
layers = view.getLayer(/ALL)
; find the savi raster and set stretch
foreach l,layers do begin
if l.data eq newRaster then l.quick_stretch='linear 5%'
endforeach
; set view to full extent
view.Zoom, /full_extent
end
Review on 12/31/2013 MM