520 Rate this article:
4.9

Automating SAR Intensity Analysis

Daryl Atencio

In my 20+ years at L3Harris Geospatial, I’ve worked with a lot of different data types and modalities, but not much with Synthetic Aperture Radar (SAR). That is until the last couple of years, when I have seen a growing interest in SAR data that I attribute in part to freely available Sentinel-1 data as well as more commercial, high-resolution SAR small satellites. So, I decided it was time to take the plunge and add SAR data to my skillset. While going through the Sentinel-1 Intensity Analysis ENVI SARscape tutorial, it occurred to me that writing a program that automates the steps would be a good learning exercise.

SAR intensity analysis is used for multiple applications such as detecting oil slicks, studying vegetation as well as monitoring vegetation over time using multiple images. Automating the extraction of the intensity product is particularly useful when extracting multiple intensity images as it increases processing speeds and reduces button clicking.

I wrote this program using IDL® 8.8.2 and ENVI® 5.6.2. You will need the “SAR Getting Started” data files, available from the ENVI Tutorials web page, to replicate this approach. To begin, extract the data to a new directory on your local drive.

To run the program:

  1. Copy the code below to a file named “sar_intensityanalysis.pro”
  2. Execute the following command at the IDL command prompt:
    IDL> sar_intensityanalysis

 

Once you run the program, you will be prompted to select the data directory. There are two input keywords to sar_intensityanalysis. See the routine’s header for information regarding these keywords.

 



;------------------------------------------------------------------------------
;+
; This helper routine allows the user to toggle between two ENVIRasterLayers
; in the ENVI display.
;-
pro SAR_ToggleLayers, eLayer, $
  DESCRIPTION=desc, $
  DESTROY=destroy, $
  END_DESCRIPTION=descEnd, $
  TITLE=title, $
  TOP_LEVEL_BASE=tlb
  compile_opt idl2, logical_predicate
  if (n_elements(eLayer) NE 2) then begin
    print, 'Invalid input: ELAYER must contain two ENVIRASTERLAYERs'
    return
  endif
  for i = 0, 1 do eLayer[i]->SetProperty, HIDE=i
  if ((n_elements(tlb) EQ 0) || ~widget_info(tlb, /VALID_ID)) then begin
    e = envi(/CURRENT)
    e->GetProperty, WIDGET_ID=wENVI
    tlb = widget_base(/COLUMN, GROUP_LEADER=wENVI, /FLOATING, $
      TITLE='Intensity Data', TLB_FRAME_ATTR=9)
    wBase = widget_base(tlb, /ALIGN_CENTER, /ROW)
    wLabel = widget_label(wBase, /DYNAMIC_RESIZE, UNAME='label', VALUE=' ')
    wBaseControl = widget_base(tlb, /COLUMN, UNAME='controls', XPAD=0, YPAD=0)
    wBase = widget_base(wBaseControl, /EXCLUSIVE, /ROW)
    wButton = [ $
      widget_button(wBase, VALUE='cross-polarized (VH)', UNAME='button_0', UVALUE=eLayer[0]), $
      widget_button(wBase, VALUE='co-polarized (VV)', UNAME='button_1', UVALUE=eLayer[1])]
    widget_control, wButton[0], /SET_BUTTON
    wBase = widget_base(wBaseControl, /ALIGN_RIGHT, /ROW)
    wButtonClose = widget_button(wBase, UNAME='button_continue', VALUE='Continue')
    widget_control, tlb, /REALIZE
  endif else begin
    widget_control, tlb, /MAP
    wBaseControl = widget_info(tlb, FIND_BY_UNAME='controls')
    wLabel = widget_info(tlb, FIND_BY_UNAME='label')
    wButton = [widget_info(tlb, FIND_BY_UNAME='button_0'), $
      widget_info(tlb, FIND_BY_UNAME='button_1')]
    wButtonClose = widget_info(tlb, FIND_BY_UNAME='button_continue')
  endelse
  widget_control, wBaseControl, /SENSITIVE
  for i = 0, 1 do widget_control, wButton[i], SET_BUTTON=(1-i)
  if (n_elements(title) EQ 1) then widget_control, tlb, BASE_SET_TITLE=title
  widget_control, wLabel, SET_VALUE=((n_elements(desc) GT 0) ? desc[0] : ' ')
  if keyword_set(destroy) then widget_control, wButtonClose, SET_VALUE='Close'
  active = !true
  while (active) do begin
    sEvent = widget_event(tlb, /NOWAIT)
    switch sEvent.id of
      wButton[0]: index = 0
      wButton[1]: begin
        if (sEvent.select EQ 1) then begin
          if (n_elements(index) EQ 0) then index = 1
          eLayer[index]->SetProperty, HIDE=0
          eLayer[1-temporary(index)]->SetProperty, HIDE=1
        endif else [] = temporary(index)
        break
      end
      
      wButtonClose: begin
        active = !false
        break
      end
      
      else:
    endswitch
  endwhile
  if keyword_set(destroy) then begin
    widget_control, tlb, /DESTROY
  endif else begin
    widget_control, wBaseControl, SENSITIVE=0
    if (n_elements(descEnd) GT 0) then begin
      widget_control, wLabel, SET_VALUE=descEnd[0]
    endif
  endelse
end

;------------------------------------------------------------------------------
;+
; This program performs all of the steps done in the Sentinel-1 Intensity
; Analysis in ENVI SARscape Tutorial
;
; :Keywords:
;   DIRECTORY: in, optional, type="string"
;     Set this to the full path of the directory containing the tutorial data.
;     If this is not set, a dialog will be presented for the user to select 
;     this folder.
;   SHOW_TOGGLE_DIALOG: in, optional, type="boolean"
;     Set this keyword to have a dialog displayed that will cause the halt 
;     execution at the end of each step.  The dialog will allow the user to 
;     toggle between the VH and VV data in the ENVI display.  
;     
;-
pro SAR_IntensityAnalysis, $
  DIRECTORY=dirRoot, $
  SHOW_TOGGLE_DIALOG=showDialog
  compile_opt idl2, logical_predicate
  err = 0
  catch, err
  if (err NE 0) then begin
    catch, /CANCEL
    help, /LAST_MESSAGE
    return
  endif
  if (n_elements(dirRoot) EQ 0) then begin
    dirRoot = dialog_pickfile(/DIRECTORY, /MUST_EXIST, $
      TITLE='Select the Tutorial Data Directory')
  endif
;
; Open and display the raw intensity data.
;
  file = filepath('sentinel1_83_20171212_095750732_IW_SIW1_D_'+['VH','VV']+'_slc_list_pwr', $
    ROOT_DIR=dirRoot, SUBDIR='2017-12-12-SLC-Power')
  if ~(total(file_test(file, /REGULAR)) EQ 2) then begin
    print, 'Input directory is missing data files'
    return
  endif
  e = envi(/CURRENT)
  if ~obj_valid(e) then e = envi()
  eView = e->GetView()
  eRaster = [e-\>OpenRaster(file[0]), e->OpenRaster(file[1])]
  eLayer = [eView->CreateLayer(eRaster[0]), eView->CreateLayer(eRaster[1])]
  eView->Zoom, LAYER=eLayer[0], /FULL_EXTENT
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayer, END_DESCRIPTION='Creating subset...', $
      TITLE='Raw Intensity Data', TOP_LEVEL_BASE=tlb
  endif
;
; Create spatial subsets of the VH and VV intensity images.
;
  fileAOI = filepath('StudySite.shp', ROOT_DIR=dirRoot)
  fileCut = filepath(file_basename(file)+'_cut', ROOT_DIR=dirRoot)  
  eTaskCut = envitask('SARsToolsManualSelection')
  eTaskCut.input_sarscapedata = file
  eTaskCut.cut_shape_file = fileAOI
  eTaskCut.output_sarscapedata = fileCut
  eTaskCut.cut_is_georeferenced = !false
  eTaskCut.USE_only_shape_corners = !true
  eTaskCut.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  eTaskCut->Execute
  eRasterCut = [e->OpenRaster(fileCut[0]), e->OpenRaster(fileCut[1])]
  eLayerCut = [eView->CreateLayer(eRasterCut[0]), $
    eView->CreateLayer(eRasterCut[1])]
  eView->Zoom, LAYER=eLayerCut[0], /FULL_EXTENT
  for i = 0, 1 do eLayer[i]->Close
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayerCut, END_DESCRIPTION='Filtering data...', $
      TITLE='Cut Intensity Data', TOP_LEVEL_BASE=tlb
  endif
;
; Apply a filter to reduce speckling noise in the VH and VV images.
;
  fileFilter = filepath(file_basename(fileCut)+'_fil', ROOT_DIR=dirRoot)
  eTaskFilter = envitask('SARsWF_ToolsGenericFilterSingleImage')
  eTaskFilter.filt_type = 'Gamma APM'
  eTaskFilter.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  for i = 0, 1 do begin
    eTaskFilter.input_sarscapedata = fileCut[i]
    eTaskFilter.output_sarscapedata = fileFilter[i]
    eTaskFilter->Execute
  endfor
  eRasterFilter = [e->OpenRaster(fileFilter[0]), e->OpenRaster(fileFilter[1])]
  eLayerFilter = [eView->CreateLayer(eRasterFilter[0]), $
    eView->CreateLayer(eRasterFilter[1])]
  for i = 0, 1 do eLayerCut[i]->SetProperty, HIDE=1
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayerFilter, $
      END_DESCRIPTION='Geocoding and radiometric correction...', $
      TITLE='Filtered Cut Data', TOP_LEVEL_BASE=tlb
  endif
;
; Create a DEM file for georeferencing.
;
  fileDEM = filepath(file_basename(fileFilter[0])+'_srtm3_dem', ROOT_DIR=dirRoot)
  eTaskDEM = envitask('SARsToolsDEMExtractionSRTM4')
  eTaskDEM.output_cartographic_system = ['GEO-GLOBAL','','GEO','','WGS84','','0.00']
  eTaskDEM.reference_sarscapedata = fileFilter
  eTaskDEM.output_sarscapedata = fileDEM
  eTaskDEM.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  eTaskDEM->\Execute
;
; Georeference the VH and VV images to a standard map projection, also called
; geocoding. As part of this step, radiometric calibration and normalization
; also will be done.
;
  fileGeo = filepath(file_basename(fileFilter)+'_geo', ROOT_DIR=dirRoot)
  eTaskGeo = envitask('SARsBasicGeocoding')
  eTaskGeo.input_sarscapedata = fileFilter
  eTaskGeo.output_sarscapedata = fileGeo
  eTaskGeo.sarscape_preference = 'Sentinel TOPSAR (IW - EW)'
  eTaskGeo.dem_sarscapedata = fileDEM
  eTaskGeo->Execute
  eRasterGeo = [e->OpenRaster(fileGeo[0]), e->OpenRaster(fileGeo[1])]
  eLayerGeo = [eView->CreateLayer(eRasterGeo[0]), $
    eView->CreateLayer(eRasterGeo[1])]
  eView->Zoom, LAYER=eLayerGeo[0], /FULL_EXTENT
  foreach eLayerClose, [eLayerCut,eLayerFilter] do eLayerClose->Close
  if keyword_set(showDialog) then begin
    SAR_ToggleLayers, eLayerGeo, /DESTROY, $
      TITLE='Georeferenced Images', TOP_LEVEL_BASE=tlb
  endif
;
; Create a color composite image from the VH and VV images.
;
  fileColor = filepath(file_basename(fileGeo[0])+'_color', ROOT_DIR=dirRoot)
  vh = eRasterGeo[0]->GetData()
  vv = eRasterGeo[1]->GetData()
  dim = size(vh, /DIMENSIONS)
  rgb = fltarr(dim[0],dim[1],3)
  rgb[*,*,0] = vv
  rgb[*,*,1] = vh
  rgb[*,*,2] = vv/vh
  eRasterColor = enviraster(rgb, URI=fileColor)
  eRasterColor->Save
  eLayerColor = eView->CreateLayer(eRasterColor)
  eView->Zoom, LAYER=eLayerColor, /FULL_EXTENT
  for i = 0, 1 do eLayerGeo[i]->SetProperty, HIDE=1
end

 

This example demonstrates calling the ENVI SARscape tasks programmatically to run the intensity analysis automatically. If you are looking for automated workflows for SAR data, ENVI SARscape Analytics walk users through some of the most common processing tasks performed using SAR data. And, if you have a problem that you think can be solved using SAR data, reach out to me or someone in our Solutions Delivery team -- we can work with you to create custom workflows based on your specific project requirements.

 

Daryl Atencio is a Solutions Delivery Engineer at L3Harris Geospatial Solutions. He received a Master of Science degree in Applied Mathematics from the University of Colorado. He started working at L3Harris in June of 2000 as a Technical Support Engineer. He joined Solutions Delivery in 2008. While in Solutions Delivery he has created custom solutions using IDL and ENVI in a wide variety of fields including non-destructive testing, mineral analysis, weather, agriculture, and tsunami forecasting.