7332
An example of performing band math with the new ENVI API
Performing band math with the new ENVI API can be a little daunting compared to the simplistic approach of the ENVI Classic routines such as MATH_DOIT. While there is some overhead associated with how to use IDL with the new ENVI API, this example makes use of the tile iterator which can be used to greatly speed up processing time and allow for processing of very large images (larger than the available RAM).
The sample program is shown below, but first there are some important bits of code that should be discussed. The first bit of code starts at the line e = ENVI() and goes to endforeach. This bit of code clears the data manager of all rasters and closes all views in ENVI. This is necessary because, when we make a new ENVI image, if the image is open in ENVI then we cannot replace the existing file. Once the data manager is clear then we can delete the file if it exists and generate a new one in it's place. This is done with the line if (file_which(thisdir, outfile) ne '') then FILE_DELETE, outfile_path.
The next important bit of code to discuss is the tile iterator itself. This process is a little different because our output raster will only have one band, so we need to specify that when we create the tiles for the tile iterator. This is how that is done:
all_tiles = (ENVIsubsetRaster(input_raster,$
BANDS= [0])).CreateTileIterator(MODE='spectral')
We need to create a temporary subset raster with only 1 band which matches our output_raster definition where NBANDS = 1. If these do not match, then the tile iterator will not work. If you need multiple output bands (say three), you can replace BANDS= [0] with BANDS= [0,1,2] as long as your original image has three bands. Because we create our tile iterator with only one band, we need to extract all of the bands from the same tile in the original image. To do this, we use the line:
all_bands = input_raster.GetData(SUB_RECT=all_tiles.CURRENT_SUBRECT)
From here, all_bands contains a 3D array of data which is an array of dimensions [1024, 1, 4]. The first dimension is set by the tile iterator, the second dimension is not used, and the third dimension contains all of the information for each of our bands. From here you can do you band math with expressions like:
b1 = all_bands[*,0,0] ;band 1
b2 = all_bands[*,0,1] ;band 2
data = (float(b1) - b2)/(b1 + b2)
Or you can directly place the values of all_bands such as:
data = (float(all_bands[*,0,0]) - all_bands[*,0,1])/(all_bands[*,0,0]+ all_bands[*,0,1])
From here, the only important thing to note is that the end of the program creates the views in ENVI and links them together. For another example of band math see the following blog post:
Developing Custom Processing Algorithms with ENVI + IDL - It’s So Easy!
To run the program, save the code below as envi_new_api_bandmath_example.pro and then you can compile and run from the IDL Workbench.
pro envi_new_api_bandmath_example
;Set compile options
compile_opt IDL2
;find the directory that this .pro file lives in so that
;we can save the output raster there
filename = ((SCOPE_TRACEBACK(/STRUCT))[-1]).FILENAME
thisdir = FILE_DIRNAME(filename)
;open ENVI, clear the data manger and views
;need to clear them because otherwise the files become locked
;we are going to replace ENVI files, so they cannot be open already
e = ENVI()
;close all open rasters
opendata = (e.data).Get(/raster,/vector)
foreach raster, opendata do raster.close
;reset our views
views = e.getview(/all)
foreach view,views do begin
;close all layers
layers = view.getlayer(/all)
foreach layer, layers do layer.close
;close all views
view.close
endforeach
;input raster that we will perform bandmath on
infile = 'qb_boulder_msi'
infile_path = FILEPATH(infile, ROOT_DIR=e.ROOT_DIR, $
SUBDIRECTORY = ['data'])
;output raster filename that we will save to
;NDI = Normalized Difference Index
outfile = 'NDI_' + infile
outfile_path = thisdir + path_sep() + outfile
;we cannot write to the output file if it exists, so we need to delete it
if (file_which(thisdir, outfile) ne '') then FILE_DELETE, outfile_path
;Load the input raster
input_raster = e.OpenRaster(infile_path)
;Set up output raster to be the same as the input input_raster
;except for number of bands. Specifically, we need
;the spatial reference to be the same
output_raster = ENVIRaster(URI = outfile_path, $
NROWS=input_raster.NROWS, $
NCOLUMNS=input_raster.NCOLUMNS, $
NBANDS=1, $
SPATIALREF=input_raster.spatialref, $
DATA_TYPE=4)
;make a tile iterator of first band for input raster because our output
;raster will only have one band
all_tiles = (ENVIsubsetRaster(input_raster,$
BANDS= [0])).CreateTileIterator(MODE='spectral')
; Iterate through tiles
FOREACH tile, all_tiles DO BEGIN
;get all spectral data from original input_raster
all_bands = input_raster.GetData(SUB_RECT=all_tiles.CURRENT_SUBRECT)
;choose our bands
b1 = all_bands[*,0,0] ;band 1
b2 = all_bands[*,0,1] ;band 2
;do the band math
data = (float(b1) - b2)/(b1 + b2)
;save the tile to the output raster
output_raster.SetTile, data, all_tiles
ENDFOREACH
;save the output_raster
output_raster.save
;print some information about where the file is saved to
print, ''
print, 'Output raster saved to:'
print, outfile_path
print,''
;load both rasters into ENVI
;dont update the ENVI display while we make our views
e.refresh, /disable
;view the original data
view1 = e.GetView()
layer1 = view1.CreateLayer(input_raster)
;view the bandmath data
view2 = e.CreateView()
layer2 = view2.CreateLayer(output_raster)
;link the two views together, make them both change when zooming
view2.GeoLink, /LINK_ALL, /ZOOM_LINK
;update ENVI display now that the views are generated and linked
e.refresh
end
reviewed by ZN 7/7/2015