3693 Rate this article:
4.2

The Shockwave Heard Around the World

Jim Pendleton

GOES-16, GOES-17, and Meteosat-8 Image Animation of the Hunga Tonga-Hunga Ha’apai Eruption’s Atmospheric Wave

Following the eruption of the Hunga Tonga-Hunga Ha’apai volcano in the South Pacific on January 15, 2022, I was inspired by some of the data products and visualizations that were initially published online. In particular, radiance difference animations of NOAA’s GOES-17 Advanced Baseline Imager (ABI, developed by L3Harris) data presented by Matthew Barlow of the University of Massachusetts Lowell on the EOS Science News website showed a very interesting phenomenon. In the animation, the Earth’s atmosphere rang like a bell due to the eruption’s shockwave, analogous to post-earthquake “ringing” of the solid Earth itself.

Reports from researchers observing ground-level barometric pressure sensors indicated that the atmospheric disturbance was detectable multiple times as the wave traversed the globe, over the course of multiple days. I wanted to extend the initial mid-level tropospheric animation in both time and spatial extent to see if that lingering effect was also visible higher in the atmosphere. (Hint: it was!)

For this example, I downloaded the “Full Disk ABI L2 Cloud and Moisture Imagery” products, limiting the imagery to band 9, only. This water vapor band is centered at 6.9 microns in the infrared and is described as the “Mid-Level Troposphere WV Band”. No surface features are visible in this wavelength band, so it was apparent continental outlines would be required for viewer orientation in the final animation.

The GOES-16 and GOES-17 satellites collect multispectral full-disk ABI image data from geostationary orbits at 75.2 and 137.2 degrees West, respectively. (An example of the GOES-17 data is shown above.)

My L3Harris Geospatial colleague here in the U.S., David Starbuck, recently published a blog that included a utility to download GOES-16 and GOES-17 data via IDL and AWS, so I will leave that exercise for the reader.

Meanwhile, the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) on EUMETSAT’s Meteosat-8 satellite also collects full-disk imagery. At the time of this event, it was in geostationary orbit at 41.5 degrees East, over the Indian Ocean. Its perspective provides coverage of an area of the Earth that is minimally visible to GOES-16 and GOES-17. Band 6 of this detector has a wavelength band centered at 7.35 microns, with a band pass that overlaps that of the GOES ABI instruments.

Because the antipode of the Hunga Tonga volcano (20.545 S, 175.394 W) is located approximately over the shared borders of Mali, Algeria, and Niger, Meteosat-8’s spatial coverage provides the opportunity to view the interference pattern of the pressure wave at the midpoint of its first traversal of the globe.

The EUMETSAT website provides access to the Meteosat-8 datasets. The option to download historical data requires registration, but you may request access via a free online account. There are multiple mechanisms provided to download data once you are logged into your account. I chose the API endpoint and the product “EO:EUM:DAT:MSG:MSG15-RSS - Rapid Scan High Rate SEVIRI Level 1.5 Image Data – MSG”, then browsed by date. A typical image file download webpage looks like the image to the left.

The image files are in EUMETSAT’s https://api.eumetsat.int/data/browse/collections “”MSG” format, with a “.nat” extension, as shown above. (Support for Meteosat MSG format file support was recently introduced in ENVI® 5.6.2). These files are quite large, so I wrote some ENVI code to extract only the single band of interest for subsequent analysis.


pro meteosat_extract_band, band_number = band_number
compile_opt strictarr
on_error, 2
band = band_number ne !null ? band_number : 6
e = envi(/current)
if (e eq !null) then e = envi(/headless)
dir = dialog_pickfile(Title = 'Select directory of Meteosat .nat files', $
    /dir, /must_exist)
if (dir eq '') then return
files = file_search(filepath(root = dir, 'MSG1-*.nat'), count = c)
if (c eq 0) then return
copybandmetadata = [ $
    'DATA GAIN VALUES', $
    'DATA OFFSET VALUES', $
    'BAND NAMES', $
    'WAVELENGTH', $
    'FWHM']
foreach file, files, index do begin
  print, (index + 1).tostring() + '/' + c.tostring(), ' ', file
  bandfile = filepath(root = file_dirname(file), $
      file_basename(file, '.nat') + '-Band-' + band.tostring() + '.dat')
  ; don't reprocess if the band file exists
  if (~file_test(bandfile)) then begin
      r1 = e.openraster(file)
      d1 = r1[0].getdata(band = band - 1)
      newraster = ENVIRaster(d1, URI = bandfile, INHERITS_FROM = r1[0])
      metadata = r1[0].metadata
      newmetadata = newraster.metadata
      foreach m, copybandmetadata do begin
          newmetadata[m] = (metadata[m])[band - 1]
      endforeach
      newraster.save
      newraster.close
  endif
  if (r1 ne !null && obj_valid(r1[0])) then foreach r, r1 do r.close
  r1 = !null
endforeach
end
	

The pressure wave is not readily apparent in any individual image from a single time. However, this feature can be enhanced by making a difference image at two adjacent time steps. The radiance difference image is literally the radiance at one time versus the radiance at the previous time step.


pro meteosat_delta_images
compile_opt strictarr
on_error, 2
e = envi(/current)
if (e eq !null) then e = envi(/headless)
 dir = dialog_pickfile($
    Title = 'Select directory of Meteosat band files for time differencing', $
    /dir, /must_exist)
if (dir eq '') then return
files = file_search(filepath(root = dir, 'MSG1-*-Band-*.dat'), count = c)
files = files.sort() ; ensure they're in time order
if (c eq 0) then return
r0 = e.openraster(files[0])
d0 = r0[0].getdata()
gain = (r0[0].metadata)['DATA GAIN VALUES']
offset = (r0[0].metadata)['DATA OFFSET VALUES']
mask = where(d0 eq 0)
d0 = temporary(d0)*gain + offset
foreach r, r0 do r.close
for i = 1, c - 1 do begin
    print, (i + 1).tostring() + '/' + c.tostring(), ' ', files[i]
    r1 = e.openraster(files[i])
    d1 = r1[0].getdata()
    mask = where(d1 eq 0)
    d1 = temporary(d1)*gain + offset
    delta = d1 - d0
    delta[mask] = 0
    fn = filepath(root = file_dirname(files[i]), $
        file_basename(files[i], '.dat') + '_diff.dat')
    if ~file_test(fn) then begin
        nr = ENVIRaster(delta, URI = fn, INHERITS_FROM = r1[0])
        nr.save
        nr.close
    endif
    foreach r, r1 do r.close
    d0 = d1
endfor
end

	  

This was repeated for the GOES ABI datasets, as well. The difference images were visually inspected to remove any “noisy” images from the sequence. Additionally, I had downloaded image sets that overlapped in time over two or more sensors, so another step involved finding the “best” image sequence to make a coherent animation, across all three sensors. I chose a sequence of GOES-17 – GOES-16 – Meteosat-8 – GOES-17 for the final video.

There are many ways to approach the task of generating an image animation for these datasets. For this example, I went old-school IDL® to generate the image overlays and annotations, using a combination of Direct Graphics and Object Graphics. Although the original imagery can be over 5,000 pixels per side, for GOES ABI data, the animation input will be normalized to only 1,024 pixels per side for all three datasets.

The following snippet grabs the pixel locations of the continental outlines for a geostationary satellite view at 41.5 degrees East for the Meteosat-8 data, in a 1024 x 1024 buffer.


window, /free, /pixmap, xsize = 1024, ysize = 1024
map_set, sat_p=[35786023.0,0,0], /satellite, 0., 41.5, /noborder ; East
map_continents, /hi
mask = tvrd(channel = 0)
wdelete, !d.window
xy = array_indices(mask, where(mask ne 0))

	  

(For the GOES data, simply modify the longitude appropriately for each satellite’s location.)

Next, for each sensor, I created an IDL SAVE-format file containing all the frames with the gray-scaled image data along with the continental overlay. They’re stored in an ordered hash with the key set to the image’s time stamp, derived from the file’s name. (Note that the NOAA and METEOSAT file-naming conventions use slightly different formats with respect to the time stamps.)

I made two datasets, one with the absolute radiance difference value clamped at 10.0, and another data set with the value clamped at 1.0. As you will see in the animation, the value at 10.0 produces “smoother” images with less noise than the 1.0 delta data, but it does not show the pressure wave as prominently or for as long.


pro generate_meteosat_frames
compile_opt strictarr
on_error, 2
dir = dialog_pickfile($
    Title = 'Select directory of Meteosat difference images', $
    /dir, /must_exist)
if (dir eq '') then return
files = file_search(filepath(root = dir, 'MSG1-*Band-6_diff.dat'), count = c)
if (c eq 0) then return
window, /free, /pixmap, xsize = 1024, ysize = 1024
map_set, sat_p=[35786023.0,0,0], /satellite, 0., 41.5, /noborder ; East
map_continents, /hi
mask = tvrd(channel = 0)
wdelete, !d.window
xy = array_indices(mask, where(mask ne 0))
deltas = ['10.', '1.']
window, /free, xsize = 1024, ysize = 1024, /pixmap
e = envi(/current)
if (e eq !null) then e = envi(/headless)
foreach delta, deltas do begin
    frames = orderedhash()
    for i = 0, c - 1 do begin
        print, (i + 1).tostring() + '/' + c.tostring(), ' ', files[i]
        r1 = e.openraster(files[i])
        d1 = r1[0].getdata()
        t = ((file_basename(files[i], '.dat')).split('-'))[5]
        tvscl, congrid($
            reverse((d1 < float(delta)) > (-float(delta)), 2), 1024, 1024)
        plots, xy[0, *], xy[1, *], psym = 3, color = 'ff0000'x, /device
        frames['c' + (t.split('\.'))[0]] = tvrd(/true)
        foreach r, r1 do r.close
    endfor
    save, frames, /compress, file = filepath(root=file_dirname(files[0]), $
        'Meteosat-8 SEVIRI Band 6_Delta_' + delta + '_Tonga.sav')
endforeach
wdelete, !d.window
end

	  

These steps were repeated for the two GOES datasets, an exercise left for the reader. The thresholded and scaled radiance difference images, with continental overlays applied, were then sequenced by time stamp over all three datasets, adding an appropriate text annotation to each.


deltas = ['10.', '1.']
foreach delta, deltas do begin
  restore, 'GOES16_ABI_Band_9_Delta_' + delta + '_Tonga.sav'
  allframes = frames
  allsources = hash()
  foreach k, frames.keys() do begin
    allsources[k] = ‘16’
  endforeach
  restore, 'GOES17 + '_ABI_Band_9_Delta_' + delta + ‘ _Tonga.sav'
  foreach k, frames.keys() do begin
    allsources[k] = ‘17’
    AllFrames[k] = frames[k]
  endforeach

  restore, 'Meteosat-8 SEVIRI Band 6_Delta_' + delta + '_degree_Tonga.sav'
  foreach k, frames.keys() do begin
    allsources[k] = ‘M’
    AllFrames[k] = frames[k]
  endforeach
  outframe = 1

  ks = (AllFrames.keys()).sort()
  time0 = 15*24+4+10/60. ; approximate eruption time
  foreach k, ks, index do begin
    source = allsources[k]
    day = strmid(k, 6, 2)
    hr = strmid(k, 8, 2)
    mi = strmid(k, 10, 2)
    t1 = day*24+hr+mi/60.
    deltat = t1 - time0
    hh = long(deltat)
    mm = abs(round((deltat - hh)*60.))
    time = 'Time from eruption (HH:MM) = ' + $
      (deltat lt 0 ? '-' : '') + hh.tostring() + ':' + string(mm, format = '(I2.2)')
    if (index eq 0) then begin
      i = image(AllFrames[k], dim=[1024, 1024], margin = [0, 0, 0, 0])
      t0  = text(.01, .95, ["Hunga Tonga-Hunga Ha'apai eruption", $
        'Sources: NOAA-NASA GOES-16 & -17 ABI L2 Band 9', $
        'Meteosat-8 SEVIRI Band ' + SEVIBand, $
        'Radiance Difference (+/- ' + delta + ')'], $
        /normal, color = 'Yellow')
      t  = text(.01, .92, time, /normal, color = 'Yellow')
    endif else begin
      i.setdata, AllFrames[k]
      t.string = time
    endelse

    write_png, filepath('frame' + string(outframe++, format = '(I3.3)') + '.png'), $
      i.copywindow()
  endforeach
  i.close
endforeach
	  

In this example, I wrote the frames individually to PNG output files in separate directories first in order to inspect them. However, the following video creation step could be integrated into that logic.


pro tongamovie
compile_opt strictarr
on_error, 2
d1 = dialog_pickfile(/must_exist, /dir, $
    title = 'Select directory of 1-delta PNG images')
if (d1 eq '') then return
d2 = dialog_pickfile(/must_exist, /dir, $
    title = 'Select directory of 10-delta PNG images')
if (d2 eq '') then return
n = file_search(d1 + '*.png', count = c)
if (c eq 0) then return
; write the MP4 video to the same directory as where this routine lives
oVid = IDLffVideoWrite(filepath(root=routine_dir(), 'tonga.mp4'), FORMAT='mp4')
fps = 3
; put the 10 and 1 radiance difference images side-by-side
newx = Round(1920/2*.75)
stream = oVid.AddVideoStream(newx*2, newx, fps)
fin = bytarr(3, 2*newx, newx)
for i = 1, c do begin
    f = 'frame' + string(i, format = '(i3.3)') + '.png'
    i1 = congrid(read_png(filepath(root = d1, f)), 3, newx, newx)
    i2 = congrid(read_png(filepath(root = d2, f)), 3, newx, newx)
    fin[0, newx, 0] = i1
    fin[0, 0, 0] = i2
    r = ovid.put(stream, fin)
endfor
end

	  

Notice that the pressure wave is still clearly visible after a full traversal back to the origin in the animation, a feature that might be quite difficult to pick out by eye in any individual, static image. I have not followed up further to see how much longer the feature is visible in the IR. I would suspect that it would require setting the radiance clamping threshold (currently set to 1.0 in this example) to lower and lower values before the signal disappears completely into the background noise as the wave’s energy slowly dissipates.

It is safe to assume that over the coming years many dissertations for advanced degrees will be generated from this singular event.

 

 

DISCOVER MORE

IDL HELPS VISUALIZE MARSHALL FIRE

PrecisionPass™

BLOG POST

See how IDL is used to pull GOES-16 ABI data, create RGB images and create an Aerosol Detection data product.

Open Source Technologies for ENVI and IDL

PrecisionPass™

RECORDED WEBINAR

Explore novel ways to use ENVI + IDL with open source technologies

DISASTER MANAGEMENT

RECORDED WEBINAR

Effectively Use Geospatial Data in the Disaster Management Cycle