利用 ENVI 与 IDL 探测闪电
原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/getting-electric-with-envi-and-idl
18274 评价本文:
3.7
利用 ENVI 与 IDL 探测闪电
Zachary Norman 2016年7月21日，星期四
几天前，我们很幸运地在我居住地（科罗拉多州博尔德附近）的东边经历了一场相当壮观的雷暴。风暴位于足够远的东方，所以我们这边没有下雨，但这段距离也让我们得以清晰地欣赏到云层中上演的精彩闪电秀。以下是我拍摄的一些风暴照片动画：
从我所在的位置可以看到大量闪电（每秒一到两次），所以我决定拿起相机拍下一些精彩的照片。拍摄完成后，这些照片就成了我这篇博客文章的数据来源。由于闪电如此频繁，我想知道是否可以利用 ENVI 的分析功能 来处理我的图像并从中检测出闪电。
事实证明，通过使用三个任务和一个栅格序列，编写代码实现这个功能相当容易。为了找到闪电，我决定使用异常检测，它根据你要寻找的目标（类似于特征提取）能很好地发现图像中的特征。通过异常检测，我找到了局部区域中与周围环境明显不同的部分，这对本场景非常适用，因为与背后的云层相比，闪电显得更亮。找到异常后，只需将异常栅格转换为分类栅格，并根据需要进行分类后处理。我进行的唯一后处理是分类筛除（classification sieving），以剔除那些被判定为异常的单像素点。整个处理过程中，我最终使用了以下任务：
为了使结果更易于可视化，我额外增加了一个步骤：将来自两个栅格序列的图像合并到一个栅格序列中。我通过提取每个栅格（原始图像和闪电像素）的数据，并生成一个包含这两个数据集的数组的栅格来实现这一点。这允许你将原始数据和经过处理的数据并排显示，从而可以看到图像中闪电应该出现的位置。
在添加了这段额外的代码之后，剩下的就是运行它。虽然花费了一些时间，但产生了一些很酷的输出。以下是闪电探测器生成的最终产品：
从动画来看，它在寻找闪电方面做得相当不错。如果闪电背后的云层非常明亮，那么闪电并不总能被找到，但当云层较暗时，闪电探测器效果很好。
如果有人想自己尝试一下，下面是我用来生成栅格序列以创建图像的 IDL 代码。你可以用自己的数据试试，然后你只需要提前创建一个栅格序列，作为过程 find_lightning（用于寻找闪电）的 INPUT_RASTER_SERIES。我使用了三个过程来生成上述输出。第一个是我用来运行另外两个过程的过程，名为 run_find_lightning。要使用它，只需将代码中'C:\path\to\raster\series'部分修改为你自己的栅格序列路径。
pro run_find_lightning compile_opt idl2
;如果ENVI尚未启动，则启动它 e = envi(/current) if (e eq !NULL) then begin e = envi(/headless) endif ;检查我们是否已经打开了栅格序列 if ISA(series, 'ENVIRASTERSERIES') then begin ;如果已打开，则回到第一个栅格 series.first endif else begin series = ENVIRasterSeries('C:\path\to\raster\series') endelse
;在图像中寻找闪电的位置 find_lightning,$ INPUT_RASTER_SERIES = series,$ OUTPUT_RASTER_SERIES = output_raster_series
;将两个栅格序列都返回到它们的第一个栅格（如果尚未如此的话） series.first output_raster_series.first
;将我们的栅格序列图像合并，生成一个超级栅格序列！ combine_series,$ INPUT_SERIES_1 = series,$ INPUT_SERIES_2 = output_raster_series,$ OUTPUT_COMBINED_SERIES = output_combined_series
end
还需要另外两个程序来调用 run_find_lightning 过程。这是第一个，find_lightning。
;寻找闪电！
pro find_lightning, $
INPUT_RASTER_SERIES = input_raster_series,$
OUTPUT_RASTER_SERIES = output_raster_series
compile_opt idl2
;获取当前ENVI会话，如果未打开则启动UI
e = envi(/current)
if (e eq !NULL) then begin
e = envi()
endif
;创建列表来保存我们的闪电栅格文件URI
lightning_raster_uris = list()
;遍历每个栅格序列
nrasters = input_raster_series.COUNT
foreach raster, input_raster_series, count do begin
;执行异常检测
AnomTask = ENVITask('RXAnomalyDetection')
AnomTask.INPUT_RASTER = raster
AnomTask.ANOMALY_DETECTION_METHOD = 'RXD'
AnomTask.MEAN_CALCULATION_METHOD = 'local'
AnomTask.KERNEL_SIZE = 15
AnomTask.execute
;打开输出
anom_raster = e.openraster(AnomTask.OUTPUT_RASTER_URI)
;通过阈值将异常转为分类
ThreshTask = ENVITask('PercentThresholdClassification')
ThreshTask.INPUT_RASTER = anom_raster
ThreshTask.THRESHOLD_PERCENT = .05
ThreshTask.execute
;打开输出
thresh_raster = e.openraster(ThreshTask.OUTPUT_RASTER_URI)
;对结果进行筛除
SieveTask = ENVITask('ClassificationSieving')
SieveTask.INPUT_RASTER = thresh_raster
SieveTask.MINIMUM_SIZE = 40
SieveTask.Execute
;打开输出
sieve_raster = e.openraster(SieveTask.OUTPUT_RASTER_URI)
;保存结果
lightning_raster_uris.add, sieve_raster.URI
;关闭中间栅格
anom_raster.close
thresh_raster.close
sieve_raster.close
;打印一些关于已处理图像数量的信息
print, string(9b) + '已完成 ' + strtrim(count+1,2) + '/' + strtrim(nrasters,2) + ' 个栅格的闪电查找'
endforeach
;将闪电栅格URI列表转换为数组
lightning_raster_uris = lightning_raster_uris.toarray()
;创建一个栅格序列
SeriesTask = ENVITask('BuildRasterSeries')
SeriesTask.INPUT_RASTER_URI = lightning_raster_uris
SeriesTask.Execute
output_raster_series = ENVIRasterSeries(SeriesTask.OUTPUT_RASTERSERIES_URI)
end
这是最后一个过程，combine_series，它将两个栅格序列合并到一个图像中。
pro combine_series,$
INPUT_SERIES_1 = input_series_1,$
INPUT_SERIES_2 = input_series_2,$
OUTPUT_COMBINED_SERIES = output_combined_series
compile_opt idl2
e = envi(/current)
if (e eq !NULL) then begin
e = envi()
endif
;保存合并图像的URI
combined_uris = list()
;将图像合并成一个超级序列
nrasters = input_series_1.COUNT
for i=0, nrasters-1 do begin
;获取栅格
image_raster = input_series_1.RASTER
lightning_raster = input_series_2.RASTER
lightning_meta = lightning_raster.metadata.dehydrate()
lightning_colors = lightning_meta['CLASS LOOKUP']
;预分配一个字节数组来保存每个栅格的数据
if (i eq 0) then begin
dims = [image_raster.NCOLUMNS, image_raster.NROWS]
data = make_array(dims[0], 2*dims[1], 3, TYPE=1, VALUE=0)
endif
;获取图像数据
image_data = image_raster.GetData(INTERLEAVE='BSQ')
;获取闪电数据
lightning_data = lightning_raster.GetData()
;将闪电数据转换为RGB
temp = reform(lightning_colors[*, lightning_data], 3, dims[0], dims[1])
;改变交错方式
lightning_data = transpose(temp, [1, 2, 0])
;用数据填充预分配数组
data[*,0:dims[1]-1,*] = image_data
data[*,dims[1]:*,*] = lightning_data
;创建新栅格
outraster = ENVIRaster(data)
outraster.save
;保存输出栅格URI
combined_uris.add, outraster.URI
;转到下一个栅格
input_series_1.next
input_series_2.next
print, string(9b) + '已合并 ' + strtrim(i+1,2) + '/' + strtrim(nrasters,2) + ' 个总栅格'
endfor
;将列表转换为数组
combined_uris = combined_uris.toarray()
;用合并后的图像创建另一个栅格序列
SeriesTask = ENVITask('BuildRasterSeries')
SeriesTask.INPUT_RASTER_URI = combined_uris
SeriesTask.Execute
;打开输出栅格序列
output_combined_series = ENVIRasterSeries(SeriesTask.OUTPUT_RASTERSERIES_URI)
end
祝您闪电狩猎愉快！