1949 Rate this article:
5.0

Getting Electric with ENVI and IDL

Zachary Norman

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!

Please login or register to post comments.