3914 Rate this article:

Automating SAR Intensity Analysis


In my 20+ years at NV5 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, $
  DESTROY=destroy, $
  TITLE=title, $
  compile_opt idl2, logical_predicate
  if (n_elements(eLayer) NE 2) then begin
    print, 'Invalid input: ELAYER must contain two ENVIRASTERLAYERs'
  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')
  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)

      wButtonClose: begin
        active = !false

  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]

; 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, $
  compile_opt idl2, logical_predicate
  err = 0
  catch, err
  if (err NE 0) then begin
    catch, /CANCEL
    help, /LAST_MESSAGE
  if (n_elements(dirRoot) EQ 0) then begin
    dirRoot = dialog_pickfile(/DIRECTORY, /MUST_EXIST, $
      TITLE='Select the Tutorial Data Directory')
  ; 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'
  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
  ; 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)'
  eRasterCut = [e->OpenRaster(fileCut[0]), e->OpenRaster(fileCut[1])]
  eLayerCut = [eView->CreateLayer(eRasterCut[0]), $
  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
  ; 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]
  eRasterFilter = [e->OpenRaster(fileFilter[0]), e->OpenRaster(fileFilter[1])]
  eLayerFilter = [eView->CreateLayer(eRasterFilter[0]), $
  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
  ; 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)'
  ; 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
  eRasterGeo = [e->OpenRaster(fileGeo[0]), e->OpenRaster(fileGeo[1])]
  eLayerGeo = [eView->CreateLayer(eRasterGeo[0]), $
  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
  ; 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)
  eLayerColor = eView->CreateLayer(eRasterColor)
  eView->Zoom, LAYER=eLayerColor, /FULL_EXTENT
  for i = 0, 1 do eLayerGeo[i]->SetProperty, HIDE=1


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 NV5 Geospatial. He received a Master of Science degree in Applied Mathematics from the University of Colorado. He started working at NV5 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.