How to find the overlap of two images with the ENVI API

There is a task called ENVIImageIntersectionTask that takes two input rasters and outputs two rasters that contain data from each input raster subset to include only the area of overlap. The following code can be used to get the image data from the intersection (overlap) and union of two images. In this example, the two rasters are subset from one larger image at the beginning of this code. But this type of thing would usually be performed using two independent images that have some overlap.

***********************

  function gridDefinitionFromRaster, oRaster
  compile_opt idl2

  oRasterCoordSys = ENVICoordSys(COORD_SYS_STR=oRaster.SPATIALREF.COORD_SYS_STR)
  gridDefinition = ENVIGridDefinition(oRasterCoordSys, $
    PIXEL_SIZE=oRaster.SPATIALREF.PIXEL_SIZE, $
    NROWS=oRaster.NROWS, $
    NCOLUMNS=oRaster.NCOLUMNS, $
    TIE_POINT_MAP=oRaster.SPATIALREF.TIE_POINT_MAP, $
    TIE_POINT_PIXEL=oRaster.SPATIALREF.TIE_POINT_PIXEL)

  return, gridDefinition

end

pro rasterIntersectionExample

  compile_opt idl2

  e=ENVI()

  ; For demonstration purposes take an original raster and take two subset
  ; rasters of different areas
  OriginalFile =  FilePath('qb_boulder_pan',SUBDIR=['data'],ROOT_DIR=e.Root_Dir)
  OriginalRaster = e.OpenRaster(OriginalFile)
  Raster1 = ENVISubsetRaster(OriginalRaster, SUB_RECT = [100,100,200,200])
  Raster2 = ENVISubsetRaster(OriginalRaster, SUB_RECT = [150,150,250,250])

  ; Get the intersection and union
  Raster1Grid = gridDefinitionFromRaster(Raster1)
  Raster2Grid = gridDefinitionFromRaster(Raster2)
  SpatialIntersection = Raster1Grid.Intersection(Raster2Grid)
  oCoordSys = ENVICoordSys(COORD_SYS_STR=Raster1.SPATIALREF.COORD_SYS_STR)
  pixelSize = Raster1.SpatialRef.Pixel_Size

  IntersectionGrid = ENVIGridDefinition(oCoordSys, EXTENT=SpatialIntersection, $

     PIXEL_SIZE=pixelSize)

  print,'Intersection', SpatialIntersection
  SpatialUnion = Raster1Grid.Union(Raster2Grid)
  print,'Union', SpatialUnion
  UnionGrid = ENVIGridDefinition(oCoordSys, EXTENT=SpatialUnion, PIXEL_SIZE=pixelSize)

  ; Create a metaspectral raster of the intersection
  Raster1Intersection = ENVISpatialGridRaster(Raster1, GRID_DEFINITION=IntersectionGrid)
  Raster2Intersection = ENVISpatialGridRaster(Raster2, GRID_DEFINITION=IntersectionGrid)
  IntersectionRaster = ENVIMetaspectralRaster([Raster1Intersection,Raster2Intersection])
  view=e.GetView()
  !null = view.CreateLayer(IntersectionRaster,BAND=0)
  !null = view.CreateLayer(IntersectionRaster,BAND=1)

  ; Create a metaspectral raster of the union
  Raster1Union = ENVISpatialGridRaster(Raster1, GRID_DEFINITION=UnionGrid)
  Raster2Union = ENVISpatialGridRaster(Raster2, GRID_DEFINITION=UnionGrid)
  UnionRaster = ENVIMetaspectralRaster([Raster1Union,Raster2Union])
  view = e.CreateView()
  !null = view.CreateLayer(UnionRaster,BAND=0)
  !null = view.CreateLayer(UnionRaster,BAND=1)

  ; Zoom to full resolution on all views
  views = e.GetView(/ALL)
  foreach view, views do begin
    layer = view.GetLayer()
    view.Zoom,/FULL_EXTENT,LAYER=layer
  endforeach

end