Getting Electric with ENVI and IDL
A few days ago we were lucky to have a pretty awesome lightning storm East of where I live (by Boulder, CO). The storm was far enough East that we weren't getting any rain, but the distance also provided an great view to the amazing lightning show that was going on in the clouds. Here is an animation of some pictures I took of the storm:
From my vantage point, you could see tons of lighting (one or two each second), so I decided to grab my camera and snap a few great pictures. After I had my pictures taken, they then became my data for this blog post. Because there was so much lightning, I wondered if I would be able to take advantage of ENVI's analytics to go through my images and detect the lightning in them.
This turned out to be pretty easy to code up with the use of three tasks and a raster series. To find the lightning, I decided to use anomaly detection which works really well to find features in images depending on what you are looking for (similar to feature extraction). With anomaly detection I found local areas that stand out from their surroundings which was perfect for this scenario because lightning is bright compared to the clouds behind it. Once you find your anomalies you then just have to turn the anomalous raster to a classification rater and perform classification cleanup if you want to. The only cleanup I performed was classification sieving to throw out single pixels which were determined to be anomalies. With all the processing, the tasks I ended up using were:
RXAnomalyDetection
PercentThresholdClassification
ClassificationSieving
In order to make the results more friendly to visualize, I took one additional step to combine the images from two raster series into one raster series. I did this by taking the data from each raster (original images and the lightning pixels) and producing one raster which had an array containing both datasets. This allows you to visualize the original data and processed data side by side so you can see where the lightning should be in the images.
After I had that extra piece of code, it just came down to running it, which took a bit of time, but produced some cool output. Here is what the lightning-detector produced for a final product:
From the animation, it does a pretty good job of finding the lightning. If the clouds are really bright behind the lightning, then the lightning isn't always found, but the lightning detector works great when the clouds are dark.
In case anyone wants to take a stab at this on their own, below is the IDL code that I used to generate the raster series used to create the images. You can give it a whirl with your own data, then you just need to create a raster series ahead of time which is the INPUT_RASTER_SERIES for the procedure find_lightning (used to find the lightning). There are three procedures I used to generate the output above. First is the procedure that I used to run my two other procedures which is called run_find_lightning. To use this, just modify the code that says 'C:\path\to\raster\series' for the path to your raster series.
pro run_find_lightning
compile_opt idl2
;start ENVI if it hasn't opened already
e = envi(/current)
if (e eq !NULL) then begin
e = envi(/headless)
endif
;check if we have opened our raster series yet
if ISA(series, 'ENVIRASTERSERIES') then begin
;return to first raster if opened already
series.first
endif else begin
series = ENVIRasterSeries('C:\path\to\raster\series')
endelse
;find where the lightning is in the images
find_lightning,$
INPUT_RASTER_SERIES = series,$
OUTPUT_RASTER_SERIES = output_raster_series
;return both raster series to their first raster in case they arent already
series.first
output_raster_series.first
;combine our raster series iamges together to produce a mega raster series!
combine_series,$
INPUT_SERIES_1 = series,$
INPUT_SERIES_2 = output_raster_series,$
OUTPUT_COMBINED_SERIES = output_combined_series
end
There are two more programs needed to call the run_find_lightning procedure. Here is the first, find lightning.
;find lightning!
pro find_lightning, $
INPUT_RASTER_SERIES = input_raster_series,$
OUTPUT_RASTER_SERIES = output_raster_series
compile_opt idl2
;get current session of ENVI, start UI if not open
e = envi(/current)
if (e eq !NULL) then begin
e = envi()
endif
;create list to hold our lightning raster files
lightning_raster_uris = list()
;iterate over each raster series
nrasters = input_raster_series.COUNT
foreach raster, input_raster_series, count do begin
;perform anomaly detection
AnomTask = ENVITask('RXAnomalyDetection')
AnomTask.INPUT_RASTER = raster
AnomTask.ANOMALY_DETECTION_METHOD = 'RXD'
AnomTask.MEAN_CALCULATION_METHOD = 'local'
AnomTask.KERNEL_SIZE = 15
AnomTask.execute
;open output
anom_raster = e.openraster(AnomTask.OUTPUT_RASTER_URI)
;threshold anomalies to classification
ThreshTask = ENVITask('PercentThresholdClassification')
ThreshTask.INPUT_RASTER = anom_raster
ThreshTask.THRESHOLD_PERCENT = .05
ThreshTask.execute
;open output
thresh_raster = e.openraster(ThreshTask.OUTPUT_RASTER_URI)
;sieve the results
SieveTask = ENVITask('ClassificationSieving')
SieveTask.INPUT_RASTER = thresh_raster
SieveTask.MINIMUM_SIZE = 40
SieveTask.Execute
;open output
sieve_raster = e.openraster(SieveTask.OUTPUT_RASTER_URI)
;save result
lightning_raster_uris.add, sieve_raster.URI
;close intermediate rasters
anom_raster.close
thresh_raster.close
sieve_raster.close
;print some info about how many images we have processed
print, string(9b) + 'Completed lightning finder on ' + strtrim(count+1,2) + ' of ' + strtrim(nrasters,2) + ' rasters'
endforeach
;convert lightning raster uris to an array
lightning_raster_uris = lightning_raster_uris.toarray()
;create a raster series
SeriesTask = ENVITask('BuildRasterSeries')
SeriesTask.INPUT_RASTER_URI = lightning_raster_uris
SeriesTask.Execute
output_raster_series = ENVIRasterSeries(SeriesTask.OUTPUT_RASTERSERIES_URI)
end
Here is the last procedure, combine_series, which combines the two raster series together in one image.
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
;save combined image URIs
combined_uris = list()
;combine the images together into one MEGA series
nrasters = input_series_1.COUNT
for i=0, nrasters-1 do begin
;get rasters
image_raster = input_series_1.RASTER
lightning_raster = input_series_2.RASTER
lightning_meta = lightning_raster.metadata.dehydrate()
lightning_colors = lightning_meta['CLASS LOOKUP']
;pre-allocate a byte array to hold the data from each raster
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
;get image data
image_data = image_raster.GetData(INTERLEAVE='BSQ')
;get ligtning data
lightning_data = lightning_raster.GetData()
;convert lightning data to RGB
temp = reform(lightning_colors[*, lightning_data], 3, dims[0], dims[1])
;change interleave
lightning_data = transpose(temp, [1, 2, 0])
;fill preallocated array with data
data[*,0:dims[1]-1,*] = image_data
data[*,dims[1]:*,*] = lightning_data
;make a new raster
outraster = ENVIRaster(data)
outraster.save
;save output raster URI
combined_uris.add, outraster.URI
;step to next rasters
input_series_1.next
input_series_2.next
print, string(9b) + 'Combined ' + strtrim(i+1,2) + ' of ' + strtrim(nrasters,2) + ' total rasters'
endfor
;convert list to array
combined_uris = combined_uris.toarray()
;create another raster series with our combined images
SeriesTask = ENVITask('BuildRasterSeries')
SeriesTask.INPUT_RASTER_URI = combined_uris
SeriesTask.Execute
;open the output raster series
output_combined_series = ENVIRasterSeries(SeriesTask.OUTPUT_RASTERSERIES_URI)
end
Happy lightning hunting!