Using R to detect the pressure wave from the 2022 Hunga Tonga eruption in personal weather station data
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:
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.
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.