The Shockwave Heard Around the World

原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/the-shockwave-heard-around-the-world

17400 Rate this article:

4.1

The Shockwave Heard Around the World

Jim Pendleton Tuesday, May 10, 2022

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 NV5 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.