Using R to detect the pressure wave from the 2022 Hunga Tonga eruption in personal weather station data

It seems like an age ago, but in fact it was only mid-January 2022 when this happened:

The satellite imagery from the Hunga Tonga eruption is unreal. Direct your attention to the lower right. The eruption then shock wave is simply incredible. pic.twitter.com/OTLCgyEozQ

— Taylor Trogdon (@TTrogdon) January 15, 2022

Wow. Now, pause for a moment and try to recall the last time you read any news about Tonga since the event.
The eruption sent an atmospheric pressure wave, clearly visible in this imagery, around the world. Friends online reported that this was detected by their personal weather stations (PWS) which made me wonder: was the wave apparent in online weather station data and can it be visualized using R?

The answers are yes and yes again.

One excellent source of weather station data is the Weather Underground. They used to have an API which could be accessed through an R package, rwunderground. This API was retired several years ago and the package no longer works. This leaves us with two options.

• find a PWS via Wundermap, access its data page (for example ISYDNE1993) and either web-scrape or copy-paste data from tables, or
• use the newer underlying API, provided via weather dot com, to access station data

Documentation for this second API is not always easy to understand or access. However, if you can obtain an API key and figure out the required endpoint, it’s quite easy to get the data into R.

Here’s a function to fetch JSON from the API and return a data frame given a PWS ID, a date and an API key.

library(tidyverse)
library(jsonlite)
library(lubridate)
theme_set(theme_bw())

json2df <- function(id, date, key) {
  u <- paste0("https://api.weather.com/v2/pws/history/all?stationId=", id, "&format=json&units=m&date=", date, "&apiKey=", key)
  u %>% 
    fromJSON() %>% 
    .$observations %>% 
    as_tibble() %>% 
    unnest(cols = c(metric))
}

We can apply the function to get data for a set of station IDs, in this case for locations near Sydney’s Northern Beaches region, like this. It assumes that the value for the API key is stored as the environment variable WUNDERGROUNDID in, for example, a .Renviron file.

apikey <- Sys.getenv('WUNDERGROUNDID')

stations <- c("ISYDNE1993", "ISYDNEY39", "ISYDNEY623", "ISYDNEY645", "ISYDNE939", "ISYDNEY181")

wudata <- lapply(stations, function(x) json2df(x, "20220115", apikey)) %>% 
  bind_rows()

The variable wudata has 37 columns. A glimpse of just the more relevant ones.

Rows: 1,728
Columns: 37
$ stationID          <chr> "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993", "ISYDNE1993",…
$ tz                 <chr> "Australia/Sydney", "Australia/Sydney", "Australia/Sydney", "Australia/Sydney", "Australia/Sydney", "Australia/Sydney", "Australia/Sydney", "Australia/Sy…
$ obsTimeUtc         <chr> "2022-01-14T13:04:58Z", "2022-01-14T13:09:46Z", "2022-01-14T13:14:50Z", "2022-01-14T13:19:54Z", "2022-01-14T13:24:58Z", "2022-01-14T13:29:46Z", "2022-01-…
$ obsTimeLocal       <chr> "2022-01-15 00:04:58", "2022-01-15 00:09:46", "2022-01-15 00:14:50", "2022-01-15 00:19:54", "2022-01-15 00:24:58", "2022-01-15 00:29:46", "2022-01-15 00:…
$ lat                <dbl> -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, -33.597, …
$ lon                <dbl> 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, 151.322, …
$ pressureMax        <dbl> 1009.41, 1009.11, 1009.31, 1009.41, 1009.21, 1009.21, 1009.21, 1009.21, 1009.11, 1009.11, 1008.94, 1009.04, 1008.74, 1008.74, 1008.74, 1008.74, 1008.64, …
$ pressureMin        <dbl> 1008.94, 1008.94, 1008.84, 1009.11, 1009.04, 1008.84, 1008.84, 1008.94, 1008.74, 1008.94, 1008.74, 1008.64, 1008.64, 1008.43, 1008.43, 1008.33, 1008.33, …
$ pressureTrend      <dbl> -5.93, 2.24, 2.54, -2.54, -0.85, -2.69, 3.39, -2.12, 1.27, -1.34, -1.27, -2.54, 0.00, -1.35, 2.54, -1.27, -2.54, 0.00, -3.81, -1.27, 2.54, -2.69, -2.54, …
[...and more...]

Air pressure is stored as pressureMax and pressureMin. Let’s use the former for the visualization.

wudata %>% 
  mutate(dt = ymd_hms(obsTimeLocal)) %>%
  ggplot(aes(dt, pressureMax)) + 
  geom_line() + 
  facet_grid(stationID ~ ., scales = "free_y") + 
  scale_x_datetime(breaks = "3 hours", date_labels = "%H:%M") +
  labs(x = "Time",
       y = "Max Pressure (hPa)",
       title = "Air pressure at 6 personal weather stations in Northern Sydney",
       subtitle = "January 15 2022") +
  theme(strip.text = element_text(size = 6))

Well, look at that. A clear spike of around 3-4 hPa. Yes, I cherry-picked some stations with the best spikes.

The rise in pressure looks like it happens at around 6 – 6:30 PM Sydney time. Does that make sense? It’s roughly 3 500 km from Sydney to Tonga. Assuming the pressure wave travels at or near the speed of sound (343 m s^-1), it should have taken 3500 / 0.343 / 60 / 60 = about 2.8 hours to reach Sydney. The eruption occurred at around 3:15 PM Sydney time. So yes: the spikes are quite consistent with a pressure wave originating from the eruption.