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:
- Copy the code below to a file named “sar_intensityanalysis.pro”
- 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.