SAR 强度分析自动化
原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/automating-sar-intensity-analysis
6082 评价本文：
4.9
SAR 强度分析自动化
作者：Daryl Atencio 日期：2022年7月5日，星期二
在我于 NV5 Geospatial 工作的 20 多年里，我接触过很多不同类型和模态的数据，但对合成孔径雷达（SAR）数据接触不多。这种情况直到最近几年才有所改变，我注意到人们对 SAR 数据的兴趣日益增长，这部分归功于免费提供的 Sentinel-1 数据以及更多商业化的高分辨率 SAR 小卫星。因此，我决定是时候深入学习，将 SAR 数据纳入我的技能范畴了。在学习 Sentinel-1 Intensity Analysis ENVI SARscape 教程时，我想到编写一个程序来自动化这些步骤将是一个很好的学习练习。
SAR 强度分析可用于多种应用，例如检测浮油、研究植被以及利用多幅图像随时间监测植被变化。在提取多幅强度图像时，自动化强度产品的提取尤其有用，因为它能提高处理速度并减少手动操作。
我使用 IDL® 8.8.2 和 ENVI® 5.6.2 编写了这个程序。要复现此方法，您需要从 ENVI 教程 网页获取“SAR Getting Started”数据文件。首先，将数据解压缩到本地驱动器的新目录中。
运行程序的步骤：
- 将以下代码复制到名为 "sar_intensityanalysis.pro" 的文件中
- 在 IDL 命令提示符下执行以下命令：
IDL> sar_intensityanalysis
运行程序后，系统会提示您选择数据目录。sar_intensityanalysis 过程有两个输入关键字。有关这些关键字的信息，请参见程序头的说明。
;------------------------------------------------------------------------------
;+
; 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
此示例演示了如何以编程方式调用 ENVI SARscape 任务来自动运行强度分析。如果您在寻找 SAR 数据的自动化工作流程，ENVI SARscape Analytics 可以引导用户完成一些使用 SAR 数据执行的最常见的处理任务。并且，如果您有一个认为可以使用 SAR 数据解决的问题，请联系我或我们 解决方案交付 团队的成员 —— 我们可以与您合作，根据您的具体项目需求创建定制的工作流程。
Daryl Atencio 是 NV5 Geospatial 的解决方案交付工程师。他获得了科罗拉多大学的应用数学硕士学位。他于 2000 年 6 月加入 NV5，担任技术支持工程师。他于 2008 年加入解决方案交付团队。在解决方案交付期间，他利用 IDL 和 ENVI 在无损检测、矿物分析、天气、农业和海啸预报等多个领域创建了定制解决方案。
