Posit AI Weblog: AO, NAO, ENSO: A wavelet research instance

Not too long ago, we confirmed use torch for wavelet research. A member of the circle of relatives of spectral research strategies, wavelet research bears some similarity to the Fourier Grow to be, and particularly, to its fashionable two-dimensional utility, the spectrogram.

As defined in that ebook excerpt, although, there are vital variations. For the needs of the present publish, it suffices to grasp that frequency-domain patterns are came upon via having a little bit “wave” (that, in reality, can also be of any form) “slide” over the information, computing stage of fit (or mismatch) locally of each and every pattern.

With this publish, then, my purpose is two-fold.

First, to introduce torchwavelets, a tiny, but helpful package deal that automates all the crucial steps concerned. In comparison to the Fourier Grow to be and its packages, the subject of wavelets is moderately “chaotic” – which means, it enjoys a lot much less shared terminology, and far much less shared apply. Because of this, it is smart for implementations to observe established, community-embraced approaches, each time such are to be had and neatly documented. With torchwavelets, we offer an implementation of Torrence and Compo’s 1998 “Sensible Information to Wavelet Research” (Torrence and Compo (1998)), an oft-cited paper that proved influential throughout a variety of utility domain names. Code-wise, our package deal can be a port of Tom Runia’s PyTorch implementation, itself in line with a previous implementation via Aaron O’Leary.

2d, to turn a lovely use case of wavelet research in a space of serious clinical pastime and super social significance (meteorology/climatology). Being in no way knowledgeable myself, I’d hope this may well be inspiring to other folks running in those fields, in addition to to scientists and analysts in different spaces the place temporal information rise up.

Concretely, what we’ll do is take 3 other atmospheric phenomena – El Niño–Southern Oscillation (ENSO), North Atlantic Oscillation (NAO), and Arctic Oscillation (AO) – and examine them the use of wavelet research. In every case, we additionally take a look at the total frequency spectrum, given via the Discrete Fourier Grow to be (DFT), in addition to a vintage time-series decomposition into development, seasonal parts, and the rest.

3 oscillations

Via a ways the best-known – probably the most notorious, I will have to say – a few of the 3 is El Niño–Southern Oscillation (ENSO), a.ok.a. El Niño/L. a. Niña. The time period refers to a converting development of sea floor temperatures and sea-level pressures going on within the equatorial Pacific. Each El Niño and L. a. Niña can and do have catastrophic affect on other folks’s lives, maximum particularly, for other folks in creating international locations west and east of the Pacific.

El Niño happens when floor water temperatures within the jap Pacific are upper than common, and the sturdy winds that in most cases blow from east to west are strangely susceptible. From April to October, this results in sizzling, extraordinarily rainy climate stipulations alongside the coasts of northern Peru and Ecuador, regularly leading to main floods. L. a. Niña, however, reasons a drop in sea floor temperatures over Southeast Asia in addition to heavy rains over Malaysia, the Philippines, and Indonesia. Whilst those are the spaces maximum gravely impacted, adjustments in ENSO reverberate around the globe.

Much less widely known than ENSO, however extremely influential as neatly, is the North Atlantic Oscillation (NAO). It strongly impacts wintry weather climate in Europe, Greenland, and North The usa. Its two states relate to the scale of the power distinction between the Icelandic Top and the Azores Low. When the power distinction is excessive, the jet movement – the ones sturdy westerly winds that blow between North The usa and Northern Europe – is but more potent than common, resulting in heat, rainy Eu winters and calmer-than-normal stipulations in Japanese North The usa. With a lower-than-normal power distinction, then again, the American East has a tendency to incur extra heavy storms and cold-air outbreaks, whilst winters in Northern Europe are chillier and extra dry.

In any case, the Arctic Oscillation (AO) is a ring-like development of sea-level power anomalies targeted on the North Pole. (Its Southern-hemisphere equal is the Antarctic Oscillation.) AO’s affect extends past the Arctic Circle, then again; it’s indicative of whether or not and what sort of Arctic air flows down into the center latitudes. AO and NAO are strongly similar, and may designate the similar bodily phenomenon at a elementary point.

Now, let’s make those characterizations extra concrete via taking a look at precise information.

Research: ENSO

We start with the best-known of those phenomena: ENSO. Knowledge are to be had from 1854 onwards; then again, for comparison with AO, we discard all information previous to January, 1950. For research, we select NINO34_MEAN, the per month moderate sea floor temperature within the Niño 3.4 area (i.e., the world between 5° South, 5° North, 190° East, and 240° East). In any case, we convert to a tsibble, the structure anticipated via feasts::STL().

library(tidyverse)
library(tsibble)

obtain.record(
  "https://bmcnoldy.rsmas.miami.edu/tropics/oni/ONI_NINO34_1854-2022.txt",
  destfile = "ONI_NINO34_1854-2022.txt"
)

enso <- read_table("ONI_NINO34_1854-2022.txt", skip = 9) %>%
  mutate(x = yearmonth(as.Date(paste0(YEAR, "-", `MON/MMM`, "-01")))) %>%
  make a selection(x, enso = NINO34_MEAN) %>%
  filter out(x >= yearmonth("1950-01"), x <= yearmonth("2022-09")) %>%
  as_tsibble(index = x)

enso

# A tsibble: 873 x 2 [1M]
          x  enso
      <mth> <dbl>
 1 1950 Jan  24.6
 2 1950 Feb  25.1
 3 1950 Mar  25.9
 4 1950 Apr  26.3
 5 1950 Would possibly  26.2
 6 1950 Jun  26.5
 7 1950 Jul  26.3
 8 1950 Aug  25.9
 9 1950 Sep  25.7
10 1950 Oct  25.7
# … with 863 extra rows

As already introduced, we need to take a look at seasonal decomposition, as neatly. Relating to seasonal periodicity, what do we predict? Except advised differently, feasts::STL() will thankfully select a window measurement for us. Alternatively, there’ll most probably be a number of essential frequencies within the information. (No longer in need of to break the suspense, however for AO and NAO, this may unquestionably be the case!). But even so, we need to compute the Fourier Grow to be anyway, so why now not do this first?

This is the ability spectrum:

Within the underneath plot, the x axis corresponds to frequencies, expressed as “choice of instances consistent with 12 months.” We most effective show frequencies as much as and together with the Nyquist frequency, i.e., part the sampling charge, which in our case is 12 (consistent with 12 months).

num_samples <- nrow(enso)
nyquist_cutoff <- ceiling(num_samples / 2) # best discernible frequency
bins_below_nyquist <- 0:nyquist_cutoff

sampling_rate <- 12 # consistent with 12 months
frequencies_per_bin <- sampling_rate / num_samples
frequencies <- frequencies_per_bin * bins_below_nyquist

df <- information.body(f = frequencies, y = as.numeric(fft[1:(nyquist_cutoff + 1)]$abs()))
df %>% ggplot(aes(f, y)) +
  geom_line() +
  xlab("frequency (consistent with 12 months)") +
  ylab("magnitude") +
  ggtitle("Spectrum of Niño 3.4 information")

There may be one dominant frequency, comparable to about annually. From this element by myself, we’d be expecting one El Niño tournament – or equivalently, one L. a. Niña – consistent with 12 months. However let’s find essential frequencies extra exactly. With now not many different periodicities status out, we would possibly as neatly limit ourselves to a few:

most powerful <- torch_topk(fft[1:(nyquist_cutoff/2)]$abs(), 3)
most powerful

[[1]]
torch_tensor
233.9855
172.2784
142.3784
[ CPUFloatType{3} ]

[[2]]
torch_tensor
74
21
7
[ CPULongType{3} ]

What now we have listed below are the magnitudes of the dominant parts, in addition to their respective boxes within the spectrum. Let’s see which precise frequencies those correspond to:

important_freqs <- frequencies[as.numeric(strongest[[2]])]
important_freqs

[1] 1.00343643 0.27491409 0.08247423 

That’s as soon as consistent with 12 months, as soon as consistent with quarter, and as soon as each and every twelve years, roughly. Or, expressed as periodicity, in the case of months (i.e., what number of months are there in a length):

num_observations_in_season <- 12/important_freqs  
num_observations_in_season

[1] 11.95890  43.65000 145.50000  

We now cross those to feasts::STL(), to procure a five-fold decomposition into development, seasonal parts, and the rest.

library(feasts)
enso %>%
  style(STL(enso ~ season(length = 12) + season(length = 44) +
              season(length = 145))) %>%
  parts() %>%
  autoplot()

In step with Loess decomposition, there nonetheless is vital noise within the information – the rest ultimate excessive in spite of our hinting at essential seasonalities. Actually, there’s no large wonder in that: Taking a look again on the DFT output, now not most effective are there many, shut to each other, low- and lowish-frequency parts, however as well as, high-frequency parts simply gained’t stop to give a contribution. And in reality, as of nowadays, ENSO forecasting – vastly essential in the case of human affect – is eager about predicting oscillation state only a 12 months prematurely. This will likely be attention-grabbing to bear in mind for once we continue to the opposite sequence – as you’ll see, it’ll most effective worsen.

Via now, we’re neatly knowledgeable about how dominant temporal rhythms resolve, or fail to resolve, what in truth occurs in environment and ocean. However we don’t know the rest about whether or not, and the way, the ones rhythms could have various in power over the time span regarded as. That is the place wavelet research is available in.

In torchwavelets, the central operation is a decision to wavelet_transform(), to instantiate an object that looks after all required operations. One argument is needed: signal_length, the choice of information issues within the sequence. And some of the defaults we want to override: dt, the time between samples, expressed within the unit we’re running with. In our case, that’s 12 months, and, having per month samples, we want to cross a price of one/12. With all different defaults untouched, research will likely be performed the use of the Morlet wavelet (to be had choices are Mexican Hat and Paul), and the turn out to be will likely be computed within the Fourier area (the quickest means, until you’ve got a GPU).

library(torchwavelets)
enso_idx <- enso$enso %>% as.numeric() %>% torch_tensor()
dt <- 1/12
wtf <- wavelet_transform(duration(enso_idx), dt = dt)

A decision to energy() will then compute the wavelet turn out to be:

power_spectrum <- wtf$energy(enso_idx)
power_spectrum$form

[1]  71 873

The result’s two-dimensional. The second one measurement holds dimension instances, i.e., the months between January, 1950 and September, 2022. The primary measurement warrants some extra rationalization.

Specifically, now we have right here the set of scales the turn out to be has been computed for. In the event you’re aware of the Fourier Grow to be and its analogue, the spectrogram, you’ll almost definitely suppose in the case of time as opposed to frequency. With wavelets, there may be an extra parameter, the size, that determines the unfold of the research development.

Some wavelets have each a scale and a frequency, through which case those can have interaction in complicated tactics. Others are outlined such that no separate frequency seems. Within the latter case, you instantly finally end up with the time vs. scale structure we see in wavelet diagrams (scaleograms). Within the former, maximum instrument hides the complexity via merging scale and frequency into one, leaving simply scale as a user-visible parameter. In torchwavelets, too, the wavelet frequency (if existent) has been “streamlined away.” Because of this, we’ll finally end up plotting time as opposed to scale, as neatly. I’ll say extra once we in truth see any such scaleogram.

For visualisation, we transpose the information and put it right into a ggplot-friendly structure:

instances <- lubridate::12 months(enso$x) + lubridate::month(enso$x) / 12
scales <- as.numeric(wtf$scales)

df <- as_tibble(as.matrix(power_spectrum$t()), .name_repair = "common") %>%
  mutate(time = instances) %>%
  pivot_longer(!time, names_to = "scale", values_to = "energy") %>%
  mutate(scale = scales[scale %>%
    str_remove("[.]{3}") %>%
    as.numeric()])
df %>% glimpse()

Rows: 61,983
Columns: 3
$ time  <dbl> 1950.083, 1950.083, 1950.083, 1950.083, 195…
$ scale <dbl> 0.1613356, 0.1759377, 0.1918614, 0.2092263,…
$ energy <dbl> 0.03617507, 0.05985500, 0.07948010, 0.09819…

There may be one further piece of knowledge to be integrated, nonetheless: the so-called “cone of affect” (COI). Visually, it is a shading that tells us which a part of the plot displays incomplete, and thus, unreliable and to-be-disregarded, information. Specifically, the larger the size, the extra spread-out the research wavelet, and the extra incomplete the overlap on the borders of the sequence when the wavelet slides over the information. You’ll see what I imply in a 2d.

The COI will get its personal information body:

And now we’re able to create the scaleogram:

labeled_scales <- c(0.25, 0.5, 1, 2, 4, 8, 16, 32, 64)
labeled_frequencies <- spherical(as.numeric(wtf$fourier_period(labeled_scales)), 1)

ggplot(df) +
  scale_y_continuous(
    trans = scales::compose_trans(scales::log2_trans(), scales::reverse_trans()),
    breaks = c(0.25, 0.5, 1, 2, 4, 8, 16, 32, 64),
    limits = c(max(scales), min(scales)),
    make bigger = c(0, 0),
    sec.axis = dup_axis(
      labels = scales::label_number(labeled_frequencies),
      title = "Fourier length (years)"
    )
  ) +
  ylab("scale (years)") +
  scale_x_continuous(breaks = seq(1950, 2020, via = 5), make bigger = c(0, 0)) +
  xlab("12 months") +
  geom_contour_filled(aes(time, scale, z = energy), display.legend = FALSE) +
  scale_fill_viridis_d(possibility = "turbo") +
  geom_ribbon(information = coi_df, aes(x = x, ymin = y, ymax = max(scales)),
              fill = "black", alpha = 0.6) +
  theme(legend.place = "none")

What we see here’s how, in ENSO, other rhythms have prevailed over the years. As a substitute of “rhythms,” I will have mentioned “scales,” or “frequencies,” or “sessions” – all the ones translate into one some other. Since, to us people, wavelet scales don’t imply that a lot, the length (in years) is displayed on an extra y axis at the proper.

So, we see that within the eighties, an (roughly) four-year length had outstanding affect. Thereafter, but longer periodicities received in dominance. And, in line with what we predict from prior research, there’s a basso continuo of annual similarity.

Additionally, be aware how, in the beginning sight, there turns out to were a decade the place a six-year length stood out: proper in the beginning of the place (for us) dimension begins, within the fifties. Alternatively, the darkish shading – the COI – tells us that, on this area, the information isn’t to be depended on.

Summing up, the two-dimensional research effectively enhances the extra compressed characterization we were given from the DFT. Ahead of we transfer directly to the following sequence, then again, let me simply temporarily cope with one query, in the event you had been questioning (if now not, simply learn on, since I gained’t be going into main points anyway): How is that this other from a spectrogram?

In a nutshell, the spectrogram splits the information into a number of “home windows,” and computes the DFT independently on they all. To compute the scaleogram, however, the research wavelet slides steadily over the information, leading to a spectrum-equivalent for the community of every pattern within the sequence. With the spectrogram, a hard and fast window measurement implies that now not all frequencies are resolved similarly neatly: The upper frequencies seem extra regularly within the period than the decrease ones, and thus, will permit for higher solution. Wavelet research, against this, is finished on a collection of scales intentionally organized in an effort to seize a huge differ of frequencies theoretically seen in a sequence of given duration.

Research: NAO

The information record for NAO is in fixed-table structure. After conversion to a tsibble, now we have:

obtain.record(
 "https://crudata.uea.ac.united kingdom/cru/information//nao/nao.dat",
 destfile = "nao.dat"
)

# wanted for AO, as neatly
use_months <- seq.Date(
  from = as.Date("1950-01-01"),
  to = as.Date("2022-09-01"),
  via = "months"
)

nao <-
  read_table(
    "nao.dat",
    col_names = FALSE,
    na = "-99.99",
    skip = 3
  ) %>%
  make a selection(-X1, -X14) %>%
  as.matrix() %>%
  t() %>%
  as.vector() %>%
  .[1:length(use_months)] %>%
  tibble(
    x = use_months,
    nao = .
  ) %>%
  mutate(x = yearmonth(x)) %>%
  fill(nao) %>%
  as_tsibble(index = x)

nao

# A tsibble: 873 x 2 [1M]
          x   nao
      <mth> <dbl>
 1 1950 Jan -0.16
 2 1950 Feb  0.25
 3 1950 Mar -1.44
 4 1950 Apr  1.46
 5 1950 Would possibly  1.34
 6 1950 Jun -3.94
 7 1950 Jul -2.75
 8 1950 Aug -0.08
 9 1950 Sep  0.19
10 1950 Oct  0.19
# … with 863 extra rows

Like prior to, we commence with the spectrum:

fft <- torch_fft_fft(as.numeric(scale(nao$nao)))

num_samples <- nrow(nao)
nyquist_cutoff <- ceiling(num_samples / 2)
bins_below_nyquist <- 0:nyquist_cutoff

sampling_rate <- 12 
frequencies_per_bin <- sampling_rate / num_samples
frequencies <- frequencies_per_bin * bins_below_nyquist

df <- information.body(f = frequencies, y = as.numeric(fft[1:(nyquist_cutoff + 1)]$abs()))
df %>% ggplot(aes(f, y)) +
  geom_line() +
  xlab("frequency (consistent with 12 months)") +
  ylab("magnitude") +
  ggtitle("Spectrum of NAO information")

Have you ever been questioning for a tiny second whether or not this was once time-domain information – now not spectral? It does glance much more noisy than the ENSO spectrum needless to say. And in reality, with NAO, predictability is way worse – forecast lead time normally quantities to only one or two weeks.

Continuing as prior to, we select dominant seasonalities (a minimum of this nonetheless is imaginable!) to cross to feasts::STL().

most powerful <- torch_topk(fft[1:(nyquist_cutoff/2)]$abs(), 6)
most powerful

[[1]]
torch_tensor
102.7191
80.5129
76.1179
75.9949
72.9086
60.8281
[ CPUFloatType{6} ]

[[2]]
torch_tensor
147
99
146
59
33
78
[ CPULongType{6} ]

important_freqs <- frequencies[as.numeric(strongest[[2]])]
important_freqs

[1] 2.0068729 1.3470790 1.9931271 0.7972509 0.4398625 1.0584192

num_observations_in_season <- 12/important_freqs  
num_observations_in_season

[1]  5.979452  8.908163  6.020690 15.051724 27.281250 11.337662

Necessary seasonal sessions are of duration six, 9, 11, fifteen, and twenty-seven months, roughly – beautiful shut in combination certainly! No surprise that, in STL decomposition, the rest is much more vital than with ENSO:

nao %>%
  style(STL(nao ~ season(length = 6) + season(length = 9) +
              season(length = 15) + season(length = 27) +
              season(length = 12))) %>%
  parts() %>%
  autoplot()

Now, what’s going to we see in the case of temporal evolution? A lot of the code that follows is equal to for ENSO, repeated right here for the reader’s comfort:

nao_idx <- nao$nao %>% as.numeric() %>% torch_tensor()
dt <- 1/12 # identical period as for ENSO
wtf <- wavelet_transform(duration(nao_idx), dt = dt)
power_spectrum <- wtf$energy(nao_idx)

instances <- lubridate::12 months(nao$x) + lubridate::month(nao$x)/12 # additionally identical
scales <- as.numeric(wtf$scales) # will likely be identical as a result of each sequence have identical duration

df <- as_tibble(as.matrix(power_spectrum$t()), .name_repair = "common") %>%
  mutate(time = instances) %>%
  pivot_longer(!time, names_to = "scale", values_to = "energy") %>%
  mutate(scale = scales[scale %>%
    str_remove("[.]{3}") %>%
    as.numeric()])

coi <- wtf$coi(instances[1], instances[length(nao_idx)])
coi_df <- information.body(x = as.numeric(coi[[1]]), y = as.numeric(coi[[2]]))

labeled_scales <- c(0.25, 0.5, 1, 2, 4, 8, 16, 32, 64) # identical since scales are identical 
labeled_frequencies <- spherical(as.numeric(wtf$fourier_period(labeled_scales)), 1)

ggplot(df) +
  scale_y_continuous(
    trans = scales::compose_trans(scales::log2_trans(), scales::reverse_trans()),
    breaks = c(0.25, 0.5, 1, 2, 4, 8, 16, 32, 64),
    limits = c(max(scales), min(scales)),
    make bigger = c(0, 0),
    sec.axis = dup_axis(
      labels = scales::label_number(labeled_frequencies),
      title = "Fourier length (years)"
    )
  ) +
  ylab("scale (years)") +
  scale_x_continuous(breaks = seq(1950, 2020, via = 5), make bigger = c(0, 0)) +
  xlab("12 months") +
  geom_contour_filled(aes(time, scale, z = energy), display.legend = FALSE) +
  scale_fill_viridis_d(possibility = "turbo") +
  geom_ribbon(information = coi_df, aes(x = x, ymin = y, ymax = max(scales)),
              fill = "black", alpha = 0.6) +
  theme(legend.place = "none")

That, actually, is a a lot more colourful image than with ENSO! Top frequencies are provide, and regularly dominant, over the entire period of time.

Apparently, although, we see similarities to ENSO, as neatly: In each, there may be the most important development, of periodicity 4 or somewhat extra years, that exerces affect all the way through the eighties, nineties, and early two-thousands – most effective with ENSO, it displays height affect all the way through the nineties, whilst with NAO, its dominance is maximum seen within the first decade of this century. Additionally, each phenomena show off a strongly seen height, of length two years, round 1970. So, is there an in depth(-ish) connection between each oscillations? This query, in fact, is for the area professionals to reply to. No less than I discovered a contemporary learn about (Scaife et al. (2014)) that now not most effective suggests there may be, however makes use of one (ENSO, the extra predictable one) to tell forecasts of the opposite:

Earlier research have proven that the El Niño–Southern Oscillation can force interannual diversifications within the NAO [Brönnimann et al., 2007] and therefore Atlantic and Eu wintry weather local weather by means of the stratosphere [Bell et al., 2009]. […] this teleconnection to the tropical Pacific is energetic in our experiments, with forecasts initialized in El Niño/L. a. Niña stipulations in November tending to be adopted via destructive/certain NAO stipulations in wintry weather.

Can we see a an identical dating for AO, our 3rd sequence beneath investigation? We may be expecting so, since AO and NAO are carefully similar (and even, two facets of the similar coin).

Research: AO

First, the information:

obtain.record(
 "https://www.cpc.ncep.noaa.gov/merchandise/precip/CWlink/daily_ao_index/per month.ao.index.b50.present.ascii.desk",
 destfile = "ao.dat"
)

ao <-
  read_table(
    "ao.dat",
    col_names = FALSE,
    skip = 1
  ) %>%
  make a selection(-X1) %>%
  as.matrix() %>% 
  t() %>%
  as.vector() %>%
  .[1:length(use_months)] %>%
  tibble(x = use_months,
         ao = .) %>%
  mutate(x = yearmonth(x)) %>%
  fill(ao) %>%
  as_tsibble(index = x) 

ao

# A tsibble: 873 x 2 [1M]
          x     ao
      <mth>  <dbl>
 1 1950 Jan -0.06 
 2 1950 Feb  0.627
 3 1950 Mar -0.008
 4 1950 Apr  0.555
 5 1950 Would possibly  0.072
 6 1950 Jun  0.539
 7 1950 Jul -0.802
 8 1950 Aug -0.851
 9 1950 Sep  0.358
10 1950 Oct -0.379
# … with 863 extra rows

And the spectrum:

fft <- torch_fft_fft(as.numeric(scale(ao$ao)))

num_samples <- nrow(ao)
nyquist_cutoff <- ceiling(num_samples / 2)
bins_below_nyquist <- 0:nyquist_cutoff

sampling_rate <- 12 # consistent with 12 months
frequencies_per_bin <- sampling_rate / num_samples
frequencies <- frequencies_per_bin * bins_below_nyquist

df <- information.body(f = frequencies, y = as.numeric(fft[1:(nyquist_cutoff + 1)]$abs()))
df %>% ggplot(aes(f, y)) +
  geom_line() +
  xlab("frequency (consistent with 12 months)") +
  ylab("magnitude") +
  ggtitle("Spectrum of AO information")

Smartly, this spectrum seems much more random than NAO’s, in that now not even a unmarried frequency stands proud. For completeness, here’s the STL decomposition:

most powerful <- torch_topk(fft[1:(nyquist_cutoff/2)]$abs(), 5)

important_freqs <- frequencies[as.numeric(strongest[[2]])]
important_freqs
# [1] 0.01374570 0.35738832 1.77319588 1.27835052 0.06872852

num_observations_in_season <- 12/important_freqs  
num_observations_in_season
# [1] 873.000000  33.576923   6.767442   9.387097 174.600000 

ao %>%
  style(STL(ao ~ season(length = 33) + season(length = 7) +
              season(length = 9) + season(length = 174))) %>%
  parts() %>%
  autoplot()

In any case, what can the scaleogram let us know about dominant patterns?

ao_idx <- ao$ao %>% as.numeric() %>% torch_tensor()
dt <- 1/12 # identical period as for ENSO and NAO
wtf <- wavelet_transform(duration(ao_idx), dt = dt)
power_spectrum <- wtf$energy(ao_idx)

instances <- lubridate::12 months(ao$x) + lubridate::month(ao$x)/12 # additionally identical
scales <- as.numeric(wtf$scales) # will likely be identical as a result of all sequence have identical duration

df <- as_tibble(as.matrix(power_spectrum$t()), .name_repair = "common") %>%
  mutate(time = instances) %>%
  pivot_longer(!time, names_to = "scale", values_to = "energy") %>%
  mutate(scale = scales[scale %>%
    str_remove("[.]{3}") %>%
    as.numeric()])

coi <- wtf$coi(instances[1], instances[length(ao_idx)])
coi_df <- information.body(x = as.numeric(coi[[1]]), y = as.numeric(coi[[2]]))

labeled_scales <- c(0.25, 0.5, 1, 2, 4, 8, 16, 32, 64) # identical since scales are identical 
labeled_frequencies <- spherical(as.numeric(wtf$fourier_period(labeled_scales)), 1)

ggplot(df) +
  scale_y_continuous(
    trans = scales::compose_trans(scales::log2_trans(), scales::reverse_trans()),
    breaks = c(0.25, 0.5, 1, 2, 4, 8, 16, 32, 64),
    limits = c(max(scales), min(scales)),
    make bigger = c(0, 0),
    sec.axis = dup_axis(
      labels = scales::label_number(labeled_frequencies),
      title = "Fourier length (years)"
    )
  ) +
  ylab("scale (years)") +
  scale_x_continuous(breaks = seq(1950, 2020, via = 5), make bigger = c(0, 0)) +
  xlab("12 months") +
  geom_contour_filled(aes(time, scale, z = energy), display.legend = FALSE) +
  scale_fill_viridis_d(possibility = "turbo") +
  geom_ribbon(information = coi_df, aes(x = x, ymin = y, ymax = max(scales)),
              fill = "black", alpha = 0.6) +
  theme(legend.place = "none")

Having observed the total spectrum, the loss of strongly dominant patterns within the scaleogram does now not come as a large wonder. It’s tempting – for me, a minimum of – to peer a mirrored image of ENSO round 1970, the entire extra since via transitivity, AO and ENSO will have to be similar one way or the other. However right here, certified judgment in reality is reserved to the professionals.

Conclusion

Like I mentioned at first, this publish could be about inspiration, now not technical element or reportable effects. And I’m hoping that inspirational it’s been, a minimum of a little bit bit. In the event you’re experimenting with wavelets your self, or plan to – or in the event you paintings within the atmospheric sciences, and wish to supply some perception at the above information/phenomena – we’d love to listen to from you!

As at all times, thank you for studying!

Photograph via ActionVance on Unsplash