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

[This article was first published on R – What You're Doing Is Rather Desperate, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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

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.

To leave a comment for the author, please follow the link and comment on their blog: R – What You're Doing Is Rather Desperate.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)