data.table::setDTthreads(1)
# Build a wave with specific wavenumber profile
waves <- list(k = 1:10,
amplitude = rnorm(10)^2,
phase = runif(10, 0, 2*pi/(1:10)))
x <- BuildWave(seq(0, 2*pi, length.out = 60)[-1], wave = waves)
# Just fancy FFT
FitWave(x, k = 1:10)
# Extract only specific wave components
plot(FilterWave(x, 1), type = "l")
plot(FilterWave(x, 2), type = "l")
plot(FilterWave(x, 1:4), type = "l")
# Remove components from the signal
plot(FilterWave(x, -4:-1), type = "l")
# The sum of the two above is the original signal (minus floating point errors)
all.equal(x, FilterWave(x, 1:4) + FilterWave(x, -4:-1))
# The Wave envelopes shows where the signal is the most "wavy".
plot(x, type = "l", col = "grey")
lines(WaveEnvelope(x), add = TRUE)
# Examples with real data
data(geopotential)
library(data.table)
# January mean of geopotential height
jan <- geopotential[month(date) == 1, .(gh = mean(gh)), by = .(lon, lat)]
# Stationary waves for each latitude
jan.waves <- jan[, FitWave(gh, 1:4), by = .(lat)]
library(ggplot2)
ggplot(jan.waves, aes(lat, amplitude, color = factor(k))) +
geom_line()
# Build field of wavenumber 1
jan[, gh.1 := BuildWave(lon*pi/180, wave = FitWave(gh, 1)), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
geom_contour(aes(z = gh.1, color = after_stat(level))) +
coord_polar()
# Build fields of wavenumber 1 and 2
waves <- jan[, BuildWave(lon*pi/180, wave = FitWave(gh, 1:2), sum = FALSE), by = .(lat)]
waves[, lon := x*180/pi]
ggplot(waves, aes(lon, lat)) +
geom_contour(aes(z = y, color = after_stat(level))) +
facet_wrap(~k) +
coord_polar()
# Field with waves 0 to 2 filtered
jan[, gh.no12 := gh - BuildWave(lon*pi/180, wave = FitWave(gh, 0:2)), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
geom_contour(aes(z = gh.no12, color = after_stat(level))) +
coord_polar()
# Much faster
jan[, gh.no12 := FilterWave(gh, -2:0), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
geom_contour(aes(z = gh.no12, color = after_stat(level))) +
coord_polar()
# Using positive numbers returns the field
jan[, gh.only12 := FilterWave(gh, 2:1), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
geom_contour(aes(z = gh.only12, color = after_stat(level))) +
coord_polar()
# Compute the envelope of the geopotential
jan[, envelope := WaveEnvelope(gh.no12), by = .(lat)]
ggplot(jan[lat == -60], aes(lon, gh.no12)) +
geom_line() +
geom_line(aes(y = envelope), color = "red")
Run the code above in your browser using DataLab